From 61cec3ca26980564e91a717179a71920b383441a Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 3 Jan 2018 21:45:04 +0900 Subject: [PATCH] Untablify --- c/_phonopy.c | 334 ++++++++++++++--------------- c/harmonic/derivative_dynmat.c | 356 +++++++++++++++---------------- c/harmonic/dynmat.c | 238 ++++++++++----------- c/harmonic_h/derivative_dynmat.h | 28 +-- c/harmonic_h/dynmat.h | 30 +-- 5 files changed, 493 insertions(+), 493 deletions(-) diff --git a/c/_phonopy.c b/c/_phonopy.c index f3da7cf3..3347886a 100644 --- a/c/_phonopy.c +++ b/c/_phonopy.c @@ -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 int distribute_fc2(double *fc2, - PHPYCONST double lat[3][3], - PHPYCONST double (*pos)[3], - const int num_pos, - const int atom_disp, - const int map_atom_disp, - PHPYCONST double r_cart[3][3], - PHPYCONST int r[3][3], - const double t[3], - const double symprec); + PHPYCONST double lat[3][3], + PHPYCONST double (*pos)[3], + const int num_pos, + const int atom_disp, + const int map_atom_disp, + PHPYCONST double r_cart[3][3], + PHPYCONST int r[3][3], + const double t[3], + const double symprec); static void distribute_fc2_with_mappings(double (*fc2)[3][3], const int * atom_list, @@ -79,11 +79,11 @@ static void distribute_fc2_with_mappings(double (*fc2)[3][3], const int num_pos); static int compute_permutation(int * rot_atom, - PHPYCONST double lat[3][3], - PHPYCONST double (*pos)[3], - PHPYCONST double (*rot_pos)[3], - const int num_pos, - const double symprec); + PHPYCONST double lat[3][3], + PHPYCONST double (*pos)[3], + PHPYCONST double (*rot_pos)[3], + const int num_pos, + const double symprec); static void gsv_copy_smallest_vectors(double (*shortest_vectors)[27][3], 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 double get_free_energy_omega(const double temperature, - const double omega); + const double omega); static double get_entropy_omega(const double temperature, - const double omega); + const double omega); 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 int check_overlap(PHPYCONST double (*pos)[3], const int num_pos, @@ -334,14 +334,14 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args) int num_satom; if (!PyArg_ParseTuple(args, "OOOOOOOO", - &dynamical_matrix, - &force_constants, - &q_vector, - &r_vector, - &multiplicity, - &mass, - &super2prim_map, - &prim2super_map)) + &dynamical_matrix, + &force_constants, + &q_vector, + &r_vector, + &multiplicity, + &mass, + &super2prim_map, + &prim2super_map)) return NULL; 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]; get_dynamical_matrix_at_q(dm, - num_patom, - num_satom, - fc, - q, - r, - multi, - m, - s2p_map, - p2s_map, - NULL, - 1); + num_patom, + num_satom, + fc, + q, + r, + multi, + m, + s2p_map, + p2s_map, + NULL, + 1); Py_RETURN_NONE; } @@ -403,17 +403,17 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args) double *charge_sum; if (!PyArg_ParseTuple(args, "OOOOOOOOOOd", - &dynamical_matrix, - &force_constants, - &q_vector, - &r_vector, - &multiplicity, - &mass, - &super2prim_map, - &prim2super_map, - &q_cart_vector, - &born, - &factor)) + &dynamical_matrix, + &force_constants, + &q_vector, + &r_vector, + &multiplicity, + &mass, + &super2prim_map, + &prim2super_map, + &q_cart_vector, + &born, + &factor)) return NULL; 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_dynamical_matrix_at_q(dm, - num_patom, - num_satom, - fc, - q, - r, - multi, - m, - s2p_map, - p2s_map, - charge_sum, - 1); + num_patom, + num_satom, + fc, + q, + r, + multi, + m, + s2p_map, + p2s_map, + charge_sum, + 1); free(charge_sum); @@ -473,14 +473,14 @@ static PyObject * py_get_dipole_dipole(PyObject *self, PyObject *args) int num_patom, num_K; if (!PyArg_ParseTuple(args, "OOOOOOOdd", - &dd_py, + &dd_py, &K_list_py, - &q_vector_py, - &q_direction_py, - &born_py, + &q_vector_py, + &q_direction_py, + &born_py, &dielectric_py, &pos_py, - &factor, + &factor, &tolerance)) return NULL; @@ -549,19 +549,19 @@ static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args) double *q_dir; if (!PyArg_ParseTuple(args, "OOOOOOOOOdOOO", - &derivative_dynmat, - &force_constants, - &q_vector, - &lattice, /* column vectors */ - &r_vector, - &multiplicity, - &mass, - &super2prim_map, - &prim2super_map, - &nac_factor, - &born, - &dielectric, - &q_direction)) { + &derivative_dynmat, + &force_constants, + &q_vector, + &lattice, /* column vectors */ + &r_vector, + &multiplicity, + &mass, + &super2prim_map, + &prim2super_map, + &nac_factor, + &born, + &dielectric, + &q_direction)) { return NULL; } @@ -594,20 +594,20 @@ static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args) } get_derivative_dynmat_at_q(ddm, - num_patom, - num_satom, - fc, - q, - lat, - r, - multi, - m, - s2p_map, - p2s_map, - nac_factor, - z, - epsilon, - q_dir); + num_patom, + num_satom, + fc, + q, + lat, + r, + multi, + m, + s2p_map, + p2s_map, + nac_factor, + z, + epsilon, + q_dir); Py_RETURN_NONE; } @@ -635,9 +635,9 @@ static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args) if (!PyArg_ParseTuple(args, "OOOO", &thermal_props_py, - &temperatures_py, - &frequencies_py, - &weights_py)) { + &temperatures_py, + &frequencies_py, + &weights_py)) { 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, - const double omega) + const double omega) { /* temperature is defined by T (K) */ /* 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, - const double omega) + const double omega) { /* temperature is defined by T (K) */ /* 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, - const double omega) + const double omega) { /* temperature is defined by T (K) */ /* 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]; distribute_fc2(fc2, - lat, - pos, - num_pos, - atom_disp, - map_atom_disp, - r_cart, - r, - t, - symprec); + lat, + pos, + num_pos, + atom_disp, + map_atom_disp, + r_cart, + r, + t, + symprec); Py_RETURN_NONE; } @@ -1070,12 +1070,12 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args) int *bz_map; if (!PyArg_ParseTuple(args, "OiOOOO", - &relative_grid_points_py, - &grid_point, - &relative_grid_address_py, - &mesh_py, - &bz_grid_address_py, - &bz_map_py)) { + &relative_grid_points_py, + &grid_point, + &relative_grid_address_py, + &mesh_py, + &bz_grid_address_py, + &bz_map_py)) { 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); thm_get_neighboring_grid_points(relative_grid_points, - grid_point, - relative_grid_address, - num_relative_grid_address, - mesh, - bz_grid_address, - bz_map); + grid_point, + relative_grid_address, + num_relative_grid_address, + mesh, + bz_grid_address, + bz_map); Py_RETURN_NONE; } @@ -1106,8 +1106,8 @@ py_thm_relative_grid_address(PyObject *self, PyObject *args) double (*reciprocal_lattice)[3]; if (!PyArg_ParseTuple(args, "OO", - &relative_grid_address_py, - &reciprocal_lattice_py)) { + &relative_grid_address_py, + &reciprocal_lattice_py)) { return NULL; } @@ -1127,7 +1127,7 @@ py_thm_all_relative_grid_address(PyObject *self, PyObject *args) int (*relative_grid_address)[24][4][3]; if (!PyArg_ParseTuple(args, "O", - &relative_grid_address_py)) { + &relative_grid_address_py)) { return NULL; } @@ -1150,9 +1150,9 @@ py_thm_integration_weight(PyObject *self, PyObject *args) double iw; if (!PyArg_ParseTuple(args, "dOs", - &omega, - &tetrahedra_omegas_py, - &function)) { + &omega, + &tetrahedra_omegas_py, + &function)) { return NULL; } @@ -1179,10 +1179,10 @@ py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args) double (*tetrahedra_omegas)[4]; if (!PyArg_ParseTuple(args, "OOOs", - &integration_weights_py, - &omegas_py, - &tetrahedra_omegas_py, - &function)) { + &integration_weights_py, + &omegas_py, + &tetrahedra_omegas_py, + &function)) { 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); thm_get_integration_weight_at_omegas(iw, - num_omegas, - omegas, - tetrahedra_omegas, - function[0]); + num_omegas, + omegas, + tetrahedra_omegas, + function[0]); Py_RETURN_NONE; } @@ -1226,13 +1226,13 @@ static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args) int address_double[3]; if (!PyArg_ParseTuple(args, "OOOOOOO", - &freq_tetras_py, - &grid_points_py, - &mesh_py, - &grid_address_py, - &gp_ir_index_py, - &relative_grid_address_py, - &frequencies_py)) { + &freq_tetras_py, + &grid_points_py, + &mesh_py, + &grid_address_py, + &gp_ir_index_py, + &relative_grid_address_py, + &frequencies_py)) { 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) for (j = 0; j < num_band * 96; j++) { for (k = 0; k < 3; k++) { - g_addr[k] = grid_address[grid_points[i]][k] + - relative_grid_address[j % 96][k]; + g_addr[k] = grid_address[grid_points[i]][k] + + relative_grid_address[j % 96][k]; } kgd_get_grid_address_double_mesh(address_double, - g_addr, - mesh, - is_shift); + g_addr, + mesh, + is_shift); gp = kgd_get_grid_point_double_mesh(address_double, mesh); 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; if (!PyArg_ParseTuple(args, "OOOOOOOO", - &dos_py, - &mesh_py, - &freq_points_py, - &frequencies_py, + &dos_py, + &mesh_py, + &freq_points_py, + &frequencies_py, &coef_py, - &grid_address_py, - &grid_mapping_table_py, - &relative_grid_address_py)) { + &grid_address_py, + &grid_mapping_table_py, + &relative_grid_address_py)) { return NULL; } @@ -1397,15 +1397,15 @@ static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args) } static int distribute_fc2(double *fc2, - PHPYCONST double lat[3][3], - PHPYCONST double (*pos)[3], - const int num_pos, - const int atom_disp, - const int map_atom_disp, - PHPYCONST double r_cart[3][3], - PHPYCONST int r[3][3], - const double t[3], - const double symprec) + PHPYCONST double lat[3][3], + PHPYCONST double (*pos)[3], + const int num_pos, + const int atom_disp, + const int map_atom_disp, + PHPYCONST double r_cart[3][3], + PHPYCONST int r[3][3], + const double t[3], + const double symprec) { int i, j, k, l, m, address_new, address; int is_found, rot_atom; @@ -1431,13 +1431,13 @@ static int distribute_fc2(double *fc2, address_new = atom_disp * num_pos * 9 + i * 9; for (j = 0; j < 3; j++) { 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][j] * r_cart[m][k] * - fc2[address + l * 3 + m]; - } - } + for (l = 0; l < 3; l++) { + for (m = 0; m < 3; m++) { + fc2[address_new + j * 3 + k] += + r_cart[l][j] * r_cart[m][k] * + fc2[address + l * 3 + m]; + } + } } } end: diff --git a/c/harmonic/derivative_dynmat.c b/c/harmonic/derivative_dynmat.c index 3791e721..4795d223 100644 --- a/c/harmonic/derivative_dynmat.c +++ b/c/harmonic/derivative_dynmat.c @@ -37,46 +37,46 @@ #define PI 3.14159265358979323846 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); + 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); + const int cart_i, + const double q[3], + const double *born); static double get_C(const double q[3], - const double *dielectric); + const double *dielectric); static double get_dA(const int atom_i, - const int cart_i, - const int cart_j, - const double *born); + 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); + const int cart_j, + const int cart_k, + const double q[3], + const double *dielectric); void get_derivative_dynmat_at_q(double *derivative_dynmat, - 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) + 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, adrs, adrsT, is_nac; 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; 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; + 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; + fabs(q[1]) < 1e-5 && + fabs(q[2]) < 1e-5) { + is_nac = 0; } } } else { @@ -108,15 +108,15 @@ void get_derivative_dynmat_at_q(double *derivative_dynmat, 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); + dnac, + num_patom, + lattice, + mass, + q, + born, + dielectric, + q_direction, + factor); } 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]); 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 (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; - } + 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 * PI); - c = cos(phase * 2 * PI); + 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 * PI); + c = cos(phase * 2 * 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 * 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 += c; + imag_phase += s; - 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 (m = 0; m < 3; m++) { + coef[m] = 0; + for (n = 0; n < 3; n++) { + coef[m] += 2 * PI * + lattice[m * 3 + n] * r[k * num_patom * 81 + i * 81 + l * 3 + n]; + } + } - 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 (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++) { - adrs = (k * num_patom * num_patom * 18 + - (i * 3 + l) * num_patom * 6 + j * 6 + m * 2); - derivative_dynmat[adrs] += ddm_real[k][l][m]; - derivative_dynmat[adrs + 1] += ddm_imag[k][l][m]; - } - } + for (l = 0; l < 3; l++) { + for (m = 0; m < 3; m++) { + adrs = (k * num_patom * num_patom * 18 + + (i * 3 + l) * num_patom * 6 + j * 6 + m * 2); + derivative_dynmat[adrs] += ddm_real[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 (j = i; j < num_patom * 3; j++) { for (k = 0; k < num_patom * 3; k++) { - 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; - derivative_dynmat[adrs] += derivative_dynmat[adrsT]; - derivative_dynmat[adrs] /= 2; - derivative_dynmat[adrs + 1] -= derivative_dynmat[adrsT+ 1]; - derivative_dynmat[adrs + 1] /= 2; - derivative_dynmat[adrsT] = derivative_dynmat[adrs]; - derivative_dynmat[adrsT + 1] = -derivative_dynmat[adrs + 1]; + 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; + derivative_dynmat[adrs] += derivative_dynmat[adrsT]; + derivative_dynmat[adrs] /= 2; + derivative_dynmat[adrs + 1] -= derivative_dynmat[adrsT+ 1]; + derivative_dynmat[adrs + 1] /= 2; + derivative_dynmat[adrsT] = derivative_dynmat[adrs]; + 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 */ /* 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) + 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; @@ -266,14 +266,14 @@ static void get_derivative_nac(double *ddnac, 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]; + q_cart[i] += rec_lat[i * 3 + j] * q_direction[j]; } 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 */ 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; - } - } - } + 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) + const int cart_i, + const double q[3], + const double *born) { int i; double sum; @@ -317,12 +317,12 @@ static double get_A(const int atom_i, 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) + const double *dielectric) { int i, j; double sum; @@ -338,33 +338,33 @@ static double get_C(const double q[3], } static double get_dA(const int atom_i, - const int cart_i, - const int cart_j, - const double *born) + 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) + 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])); + 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])); + 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])); + q[0] * (dielectric[2] + dielectric[6]) + + q[1] * (dielectric[5] + dielectric[7])); } return 0; } diff --git a/c/harmonic/dynmat.c b/c/harmonic/dynmat.c index 133a0500..b8347f68 100644 --- a/c/harmonic/dynmat.c +++ b/c/harmonic/dynmat.c @@ -38,83 +38,83 @@ #define PI 3.14159265358979323846 static void get_dynmat_ij(double *dynamical_matrix, - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double *mass, - const int *s2p_map, - const int *p2s_map, - const double *charge_sum, - const int i, - const int j); + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double *mass, + const int *s2p_map, + const int *p2s_map, + const double *charge_sum, + const int i, + const int j); static void get_dm(double dm_real[3][3], - double dm_imag[3][3], - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double mass_sqrt, - const int *p2s_map, - const double *charge_sum, - const int i, - const int j, - const int k); + double dm_imag[3][3], + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double mass_sqrt, + const int *p2s_map, + const double *charge_sum, + const int i, + const int j, + const int k); static double get_dielectric_part(const double q[3], const double *dielectric); int get_dynamical_matrix_at_q(double *dynamical_matrix, - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double *mass, - const int *s2p_map, - const int *p2s_map, - const double *charge_sum, - const int with_openmp) + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double *mass, + const int *s2p_map, + const int *p2s_map, + const double *charge_sum, + const int with_openmp) { int i, j, ij, adrs, adrsT; - + if (with_openmp) { #pragma omp parallel for for (ij = 0; ij < num_patom * num_patom ; ij++) { get_dynmat_ij(dynamical_matrix, - num_patom, - num_satom, - fc, - q, - r, - multi, - mass, - s2p_map, - p2s_map, - charge_sum, - ij / num_patom, /* i */ - ij % num_patom); /* j */ + num_patom, + num_satom, + fc, + q, + r, + multi, + mass, + s2p_map, + p2s_map, + charge_sum, + ij / num_patom, /* i */ + ij % num_patom); /* j */ } } else { for (i = 0; i < num_patom; i++) { for (j = 0; j < num_patom; j++) { - get_dynmat_ij(dynamical_matrix, - num_patom, - num_satom, - fc, - q, - r, - multi, - mass, - s2p_map, - p2s_map, - charge_sum, - i, - j); + get_dynmat_ij(dynamical_matrix, + num_patom, + num_satom, + fc, + q, + r, + multi, + mass, + s2p_map, + p2s_map, + charge_sum, + i, + j); } } } @@ -170,10 +170,10 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */ if (!q_direction) { continue; } else { - for (i = 0; i < 3; i++) {q_K[i] = q_direction[i];} + for (i = 0; i < 3; i++) {q_K[i] = q_direction[i];} } } else { - for (i = 0; i < 3; i++) {q_K[i] = K_list[g * 3 + i];} + for (i = 0; i < 3; i++) {q_K[i] = K_list[g * 3 + i];} } get_charge_sum(charge_sum, @@ -247,14 +247,14 @@ void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */ } void get_charge_sum(double *charge_sum, - const int num_patom, - const double factor, /* 4pi/V*unit-conv and denominator */ - const double q_vector[3], - const double *born) + const int num_patom, + const double factor, /* 4pi/V*unit-conv and denominator */ + const double q_vector[3], + const double *born) { int i, j, k, a, b; double (*q_born)[3]; - + q_born = (double (*)[3]) malloc(sizeof(double[3]) * num_patom); for (i = 0; i < num_patom; i++) { for (j = 0; j < 3; j++) { @@ -265,7 +265,7 @@ void get_charge_sum(double *charge_sum, for (i = 0; i < num_patom; i++) { for (j = 0; j < 3; j++) { 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 (j = 0; j < num_patom; j++) { for (a = 0; a < 3; a++) { - for (b = 0; b < 3; b++) { - charge_sum[i * 9 * num_patom + j * 9 + a * 3 + b] = - q_born[i][a] * q_born[j][b] * factor; - } + for (b = 0; b < 3; b++) { + charge_sum[i * 9 * num_patom + j * 9 + a * 3 + b] = + q_born[i][a] * q_born[j][b] * factor; + } } } } @@ -287,23 +287,23 @@ void get_charge_sum(double *charge_sum, static void get_dynmat_ij(double *dynamical_matrix, - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double *mass, - const int *s2p_map, - const int *p2s_map, - const double *charge_sum, - const int i, - const int j) + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double *mass, + const int *s2p_map, + const int *p2s_map, + const double *charge_sum, + const int i, + const int j) { int k, l, adrs; double mass_sqrt; double dm_real[3][3], dm_imag[3][3]; - + mass_sqrt = sqrt(mass[i] * mass[j]); for (k = 0; k < 3; k++) { @@ -318,21 +318,21 @@ static void get_dynmat_ij(double *dynamical_matrix, continue; } get_dm(dm_real, - dm_imag, - num_patom, - num_satom, - fc, - q, - r, - multi, - mass_sqrt, - p2s_map, - charge_sum, - i, - j, - k); + dm_imag, + num_patom, + num_satom, + fc, + q, + r, + multi, + mass_sqrt, + p2s_map, + charge_sum, + i, + j, + k); } - + for (k = 0; k < 3; k++) { for (l = 0; l < 3; l++) { adrs = (i * 3 + k) * num_patom * 6 + j * 6 + l * 2; @@ -343,19 +343,19 @@ static void get_dynmat_ij(double *dynamical_matrix, } static void get_dm(double dm_real[3][3], - double dm_imag[3][3], - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double mass_sqrt, - const int *p2s_map, - const double *charge_sum, - const int i, - const int j, - const int k) + double dm_imag[3][3], + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double mass_sqrt, + const int *p2s_map, + const double *charge_sum, + const int i, + const int j, + const int k) { int l, m; 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 (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; + fc_elem = fc[p2s_map[i] * num_satom * 9 + + k * 9 + l * 3 + m] / mass_sqrt; } dm_real[l][m] += fc_elem * cos_phase; dm_imag[l][m] += fc_elem * sin_phase; diff --git a/c/harmonic_h/derivative_dynmat.h b/c/harmonic_h/derivative_dynmat.h index de0d1a67..00547135 100644 --- a/c/harmonic_h/derivative_dynmat.h +++ b/c/harmonic_h/derivative_dynmat.h @@ -36,19 +36,19 @@ #define __derivative_dynmat_H__ void get_derivative_dynmat_at_q(double *derivative_dynmat, - 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); + 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 diff --git a/c/harmonic_h/dynmat.h b/c/harmonic_h/dynmat.h index d90a6d71..0ca03fcc 100644 --- a/c/harmonic_h/dynmat.h +++ b/c/harmonic_h/dynmat.h @@ -36,17 +36,17 @@ #define __dynmat_H__ int get_dynamical_matrix_at_q(double *dynamical_matrix, - const int num_patom, - const int num_satom, - const double *fc, - const double *q, - const double *r, - const int *multi, - const double *mass, - const int *s2p_map, - const int *p2s_map, - const double *charge_sum, - const int with_openmp); + const int num_patom, + const int num_satom, + const double *fc, + const double *q, + const double *r, + const int *multi, + const double *mass, + const int *s2p_map, + const int *p2s_map, + const double *charge_sum, + const int with_openmp); void get_dipole_dipole(double *dd, /* [natom, 3, natom, 3, (real, imag)] */ const double *K_list, /* [num_kvec, 3] */ 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 tolerance); void get_charge_sum(double *charge_sum, - const int num_patom, - const double factor, - const double q_vector[3], - const double *born); + const int num_patom, + const double factor, + const double q_vector[3], + const double *born); #endif