1 /* Copyright (C) 2011 Atsushi Togo */
2 /* All rights reserved. */
3 
4 /* This file is part of phonopy. */
5 
6 /* Redistribution and use in source and binary forms, with or without */
7 /* modification, are permitted provided that the following conditions */
8 /* are met: */
9 
10 /* * Redistributions of source code must retain the above copyright */
11 /*   notice, this list of conditions and the following disclaimer. */
12 
13 /* * Redistributions in binary form must reproduce the above copyright */
14 /*   notice, this list of conditions and the following disclaimer in */
15 /*   the documentation and/or other materials provided with the */
16 /*   distribution. */
17 
18 /* * Neither the name of the phonopy project nor the names of its */
19 /*   contributors may be used to endorse or promote products derived */
20 /*   from this software without specific prior written permission. */
21 
22 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
23 /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
24 /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
25 /* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
26 /* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
27 /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
28 /* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
29 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
30 /* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
31 /* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
32 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
33 /* POSSIBILITY OF SUCH DAMAGE. */
34 
35 #include <Python.h>
36 #include <stdio.h>
37 #include <stddef.h>
38 #include <math.h>
39 #include <float.h>
40 #include <numpy/arrayobject.h>
41 #include "phonopy.h"
42 
43 /* PHPYCONST is defined in dynmat.h */
44 
45 /* Build dynamical matrix */
46 static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args);
47 static PyObject * py_perm_trans_symmetrize_fc(PyObject *self, PyObject *args);
48 static PyObject *
49 py_perm_trans_symmetrize_compact_fc(PyObject *self, PyObject *args);
50 static PyObject * py_transpose_compact_fc(PyObject *self, PyObject *args);
51 static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args);
52 static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args);
53 static PyObject * py_get_recip_dipole_dipole(PyObject *self, PyObject *args);
54 static PyObject * py_get_recip_dipole_dipole_q0(PyObject *self, PyObject *args);
55 static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args);
56 static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args);
57 static PyObject * py_distribute_fc2(PyObject *self, PyObject *args);
58 static PyObject * py_compute_permutation(PyObject *self, PyObject *args);
59 static PyObject *
60 py_gsv_set_smallest_vectors_sparse(PyObject *self, PyObject *args);
61 static PyObject *
62 py_gsv_set_smallest_vectors_dense(PyObject *self, PyObject *args);
63 static PyObject *
64 py_thm_relative_grid_address(PyObject *self, PyObject *args);
65 static PyObject *
66 py_thm_all_relative_grid_address(PyObject *self, PyObject *args);
67 static PyObject *
68 py_thm_integration_weight(PyObject *self, PyObject *args);
69 static PyObject *
70 py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args);
71 static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args);
72 static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args);
73 
74 struct module_state {
75   PyObject *error;
76 };
77 
78 #if PY_MAJOR_VERSION >= 3
79 #define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))
80 #else
81 #define GETSTATE(m) (&_state)
82 static struct module_state _state;
83 #endif
84 
85 static PyObject *
error_out(PyObject * m)86 error_out(PyObject *m) {
87   struct module_state *st = GETSTATE(m);
88   PyErr_SetString(st->error, "something bad happened");
89   return NULL;
90 }
91 
92 static PyMethodDef _phonopy_methods[] = {
93   {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
94   {"transform_dynmat_to_fc", py_transform_dynmat_to_fc, METH_VARARGS,
95    "Transform a set of dynmat to force constants"},
96   {"perm_trans_symmetrize_fc", py_perm_trans_symmetrize_fc, METH_VARARGS,
97    "Enforce permutation and translational symmetry of force constants"},
98   {"perm_trans_symmetrize_compact_fc", py_perm_trans_symmetrize_compact_fc,
99    METH_VARARGS,
100    "Enforce permutation and translational symmetry of compact force constants"},
101   {"transpose_compact_fc", py_transpose_compact_fc,
102    METH_VARARGS,
103    "Transpose compact force constants"},
104   {"dynamical_matrix", py_get_dynamical_matrix, METH_VARARGS,
105    "Dynamical matrix"},
106   {"nac_dynamical_matrix", py_get_nac_dynamical_matrix, METH_VARARGS,
107    "NAC dynamical matrix"},
108   {"recip_dipole_dipole", py_get_recip_dipole_dipole, METH_VARARGS,
109    "Reciprocal part of dipole-dipole interaction"},
110   {"recip_dipole_dipole_q0", py_get_recip_dipole_dipole_q0, METH_VARARGS,
111    "q=0 terms of reciprocal part of dipole-dipole interaction"},
112   {"derivative_dynmat", py_get_derivative_dynmat, METH_VARARGS,
113    "Q derivative of dynamical matrix"},
114   {"thermal_properties", py_get_thermal_properties, METH_VARARGS,
115    "Thermal properties"},
116   {"distribute_fc2", py_distribute_fc2, METH_VARARGS,
117    "Distribute force constants for all atoms in atom_list using precomputed symmetry mappings."},
118   {"compute_permutation", py_compute_permutation, METH_VARARGS,
119    "Compute indices of original points in a set of rotated points."},
120   {"gsv_set_smallest_vectors_sparse", py_gsv_set_smallest_vectors_sparse,
121    METH_VARARGS, "Set shortest vectors in sparse array."},
122   {"gsv_set_smallest_vectors_dense", py_gsv_set_smallest_vectors_dense,
123    METH_VARARGS, "Set shortest vectors in dense array."},
124   {"tetrahedra_relative_grid_address", py_thm_relative_grid_address,
125    METH_VARARGS, "Relative grid addresses of vertices of 24 tetrahedra"},
126   {"all_tetrahedra_relative_grid_address",
127    py_thm_all_relative_grid_address, METH_VARARGS,
128    "4 (all) sets of relative grid addresses of vertices of 24 tetrahedra"},
129   {"tetrahedra_integration_weight", py_thm_integration_weight,
130    METH_VARARGS, "Integration weight for tetrahedron method"},
131   {"tetrahedra_integration_weight_at_omegas",
132    py_thm_integration_weight_at_omegas,
133    METH_VARARGS, "Integration weight for tetrahedron method at omegas"},
134   {"tetrahedra_frequencies", py_get_tetrahedra_frequenies,
135    METH_VARARGS, "Run tetrahedron method"},
136   {"tetrahedron_method_dos", py_tetrahedron_method_dos,
137    METH_VARARGS, "Run tetrahedron method"},
138   {NULL, NULL, 0, NULL}
139 };
140 
141 #if PY_MAJOR_VERSION >= 3
142 
_phonopy_traverse(PyObject * m,visitproc visit,void * arg)143 static int _phonopy_traverse(PyObject *m, visitproc visit, void *arg) {
144   Py_VISIT(GETSTATE(m)->error);
145   return 0;
146 }
147 
_phonopy_clear(PyObject * m)148 static int _phonopy_clear(PyObject *m) {
149   Py_CLEAR(GETSTATE(m)->error);
150   return 0;
151 }
152 
153 static struct PyModuleDef moduledef = {
154   PyModuleDef_HEAD_INIT,
155   "_phonopy",
156   NULL,
157   sizeof(struct module_state),
158   _phonopy_methods,
159   NULL,
160   _phonopy_traverse,
161   _phonopy_clear,
162   NULL
163 };
164 
165 #define INITERROR return NULL
166 
167 PyObject *
PyInit__phonopy(void)168 PyInit__phonopy(void)
169 
170 #else
171 #define INITERROR return
172 
173   void
174   init_phonopy(void)
175 #endif
176 {
177 #if PY_MAJOR_VERSION >= 3
178   PyObject *module = PyModule_Create(&moduledef);
179 #else
180   PyObject *module = Py_InitModule("_phonopy", _phonopy_methods);
181 #endif
182   struct module_state *st;
183   if (module == NULL)
184     INITERROR;
185   st = GETSTATE(module);
186 
187   st->error = PyErr_NewException("_phonopy.Error", NULL, NULL);
188   if (st->error == NULL) {
189     Py_DECREF(module);
190     INITERROR;
191   }
192 
193 #if PY_MAJOR_VERSION >= 3
194   return module;
195 #endif
196 }
197 
py_transform_dynmat_to_fc(PyObject * self,PyObject * args)198 static PyObject * py_transform_dynmat_to_fc(PyObject *self, PyObject *args)
199 {
200   PyArrayObject* py_force_constants;
201   PyArrayObject* py_dynamical_matrices;
202   PyArrayObject* py_commensurate_points;
203   PyArrayObject* py_svecs;
204   PyArrayObject* py_multi;
205   PyArrayObject* py_masses;
206   PyArrayObject* py_s2pp_map;
207   PyArrayObject* py_fc_index_map;
208 
209   double* fc;
210   double* dm;
211   double (*comm_points)[3];
212   double (*svecs)[3];
213   double* masses;
214   long (*multi)[2];
215   long* s2pp_map;
216   long* fc_index_map;
217   long num_patom;
218   long num_satom;
219 
220   if (!PyArg_ParseTuple(args, "OOOOOOOO",
221                         &py_force_constants,
222                         &py_dynamical_matrices,
223                         &py_commensurate_points,
224                         &py_svecs,
225                         &py_multi,
226                         &py_masses,
227                         &py_s2pp_map,
228                         &py_fc_index_map)) {
229     return NULL;
230   }
231 
232   fc = (double*)PyArray_DATA(py_force_constants);
233   dm = (double*)PyArray_DATA(py_dynamical_matrices);
234   comm_points = (double(*)[3])PyArray_DATA(py_commensurate_points);
235   svecs = (double(*)[3])PyArray_DATA(py_svecs);
236   masses = (double*)PyArray_DATA(py_masses);
237   multi = (long(*)[2])PyArray_DATA(py_multi);
238   s2pp_map = (long*)PyArray_DATA(py_s2pp_map);
239   fc_index_map = (long*)PyArray_DATA(py_fc_index_map);
240   num_patom = PyArray_DIMS(py_multi)[1];
241   num_satom = PyArray_DIMS(py_multi)[0];
242 
243   phpy_transform_dynmat_to_fc(fc,
244                               dm,
245                               comm_points,
246                               svecs,
247                               multi,
248                               masses,
249                               s2pp_map,
250                               fc_index_map,
251                               num_patom,
252                               num_satom);
253 
254   Py_RETURN_NONE;
255 }
256 
py_compute_permutation(PyObject * self,PyObject * args)257 static PyObject * py_compute_permutation(PyObject *self, PyObject *args)
258 {
259   PyArrayObject* permutation;
260   PyArrayObject* lattice;
261   PyArrayObject* positions;
262   PyArrayObject* permuted_positions;
263   double symprec;
264 
265   int* rot_atoms;
266   double (*lat)[3];
267   double (*pos)[3];
268   double (*rot_pos)[3];
269   int num_pos;
270 
271   int is_found;
272 
273   if (!PyArg_ParseTuple(args, "OOOOd",
274                         &permutation,
275                         &lattice,
276                         &positions,
277                         &permuted_positions,
278                         &symprec)) {
279     return NULL;
280   }
281 
282   rot_atoms = (int*)PyArray_DATA(permutation);
283   lat = (double(*)[3])PyArray_DATA(lattice);
284   pos = (double(*)[3])PyArray_DATA(positions);
285   rot_pos = (double(*)[3])PyArray_DATA(permuted_positions);
286   num_pos = PyArray_DIMS(positions)[0];
287 
288   is_found = phpy_compute_permutation(rot_atoms,
289                                       lat,
290                                       pos,
291                                       rot_pos,
292                                       num_pos,
293                                       symprec);
294 
295   if (is_found) {
296     Py_RETURN_TRUE;
297   } else {
298     Py_RETURN_FALSE;
299   }
300 }
301 
302 static PyObject *
py_gsv_set_smallest_vectors_sparse(PyObject * self,PyObject * args)303 py_gsv_set_smallest_vectors_sparse(PyObject *self, PyObject *args)
304 {
305   PyArrayObject* py_smallest_vectors;
306   PyArrayObject* py_multiplicity;
307   PyArrayObject* py_pos_to;
308   PyArrayObject* py_pos_from;
309   PyArrayObject* py_lattice_points;
310   PyArrayObject* py_reduced_basis;
311   PyArrayObject* py_trans_mat;
312   double symprec;
313 
314   double (*smallest_vectors)[27][3];
315   int * multiplicity;
316   double (*pos_to)[3];
317   double (*pos_from)[3];
318   int (*lattice_points)[3];
319   double (*reduced_basis)[3];
320   int (*trans_mat)[3];
321   int num_pos_to, num_pos_from, num_lattice_points;
322 
323   if (!PyArg_ParseTuple(args, "OOOOOOOd",
324                         &py_smallest_vectors,
325                         &py_multiplicity,
326                         &py_pos_to,
327                         &py_pos_from,
328                         &py_lattice_points,
329                         &py_reduced_basis,
330                         &py_trans_mat,
331                         &symprec)) {
332     return NULL;
333   }
334 
335   smallest_vectors = (double(*)[27][3])PyArray_DATA(py_smallest_vectors);
336   multiplicity = (int*)PyArray_DATA(py_multiplicity);
337   pos_to = (double(*)[3])PyArray_DATA(py_pos_to);
338   pos_from = (double(*)[3])PyArray_DATA(py_pos_from);
339   num_pos_to = PyArray_DIMS(py_pos_to)[0];
340   num_pos_from = PyArray_DIMS(py_pos_from)[0];
341   lattice_points = (int(*)[3])PyArray_DATA(py_lattice_points);
342   num_lattice_points = PyArray_DIMS(py_lattice_points)[0];
343   reduced_basis = (double(*)[3])PyArray_DATA(py_reduced_basis);
344   trans_mat = (int(*)[3])PyArray_DATA(py_trans_mat);
345 
346   phpy_set_smallest_vectors_sparse(smallest_vectors,
347                                    multiplicity,
348                                    pos_to,
349                                    num_pos_to,
350                                    pos_from,
351                                    num_pos_from,
352                                    lattice_points,
353                                    num_lattice_points,
354                                    reduced_basis,
355                                    trans_mat,
356                                    symprec);
357 
358   Py_RETURN_NONE;
359 }
360 
361 static PyObject *
py_gsv_set_smallest_vectors_dense(PyObject * self,PyObject * args)362 py_gsv_set_smallest_vectors_dense(PyObject *self, PyObject *args)
363 {
364   PyArrayObject* py_smallest_vectors;
365   PyArrayObject* py_multiplicity;
366   PyArrayObject* py_pos_to;
367   PyArrayObject* py_pos_from;
368   PyArrayObject* py_lattice_points;
369   PyArrayObject* py_reduced_basis;
370   PyArrayObject* py_trans_mat;
371   long initialize;
372   double symprec;
373 
374   double (*smallest_vectors)[3];
375   long (*multiplicity)[2];
376   double (*pos_to)[3];
377   double (*pos_from)[3];
378   long (*lattice_points)[3];
379   double (*reduced_basis)[3];
380   long (*trans_mat)[3];
381   long num_pos_to, num_pos_from, num_lattice_points;
382 
383   if (!PyArg_ParseTuple(args, "OOOOOOOld",
384                         &py_smallest_vectors,
385                         &py_multiplicity,
386                         &py_pos_to,
387                         &py_pos_from,
388                         &py_lattice_points,
389                         &py_reduced_basis,
390                         &py_trans_mat,
391                         &initialize,
392                         &symprec)) {
393     return NULL;
394   }
395 
396   smallest_vectors = (double(*)[3])PyArray_DATA(py_smallest_vectors);
397   multiplicity = (long(*)[2])PyArray_DATA(py_multiplicity);
398   pos_to = (double(*)[3])PyArray_DATA(py_pos_to);
399   pos_from = (double(*)[3])PyArray_DATA(py_pos_from);
400   num_pos_to = PyArray_DIMS(py_pos_to)[0];
401   num_pos_from = PyArray_DIMS(py_pos_from)[0];
402   lattice_points = (long(*)[3])PyArray_DATA(py_lattice_points);
403   num_lattice_points = PyArray_DIMS(py_lattice_points)[0];
404   reduced_basis = (double(*)[3])PyArray_DATA(py_reduced_basis);
405   trans_mat = (long(*)[3])PyArray_DATA(py_trans_mat);
406 
407   phpy_set_smallest_vectors_dense(smallest_vectors,
408                                   multiplicity,
409                                   pos_to,
410                                   num_pos_to,
411                                   pos_from,
412                                   num_pos_from,
413                                   lattice_points,
414                                   num_lattice_points,
415                                   reduced_basis,
416                                   trans_mat,
417                                   initialize,
418                                   symprec);
419 
420   Py_RETURN_NONE;
421 }
422 
py_perm_trans_symmetrize_fc(PyObject * self,PyObject * args)423 static PyObject * py_perm_trans_symmetrize_fc(PyObject *self, PyObject *args)
424 {
425   PyArrayObject* py_fc;
426   double *fc;
427   int level;
428 
429   int n_satom;
430 
431   if (!PyArg_ParseTuple(args, "Oi", &py_fc, &level)) {
432     return NULL;
433   }
434 
435   fc = (double*)PyArray_DATA(py_fc);
436   n_satom = PyArray_DIMS(py_fc)[0];
437 
438   phpy_perm_trans_symmetrize_fc(fc, n_satom, level);
439 
440   Py_RETURN_NONE;
441 }
442 
443 static PyObject *
py_perm_trans_symmetrize_compact_fc(PyObject * self,PyObject * args)444 py_perm_trans_symmetrize_compact_fc(PyObject *self, PyObject *args)
445 {
446   PyArrayObject* py_fc;
447   PyArrayObject* py_permutations;
448   PyArrayObject* py_s2pp_map;
449   PyArrayObject* py_p2s_map;
450   PyArrayObject* py_nsym_list;
451   int level;
452   double *fc;
453   int *perms;
454   int *s2pp;
455   int *p2s;
456   int *nsym_list;
457 
458   int n_patom, n_satom;
459 
460   if (!PyArg_ParseTuple(args, "OOOOOi",
461                         &py_fc,
462                         &py_permutations,
463                         &py_s2pp_map,
464                         &py_p2s_map,
465                         &py_nsym_list,
466                         &level)) {
467     return NULL;
468   }
469 
470   fc = (double*)PyArray_DATA(py_fc);
471   perms = (int*)PyArray_DATA(py_permutations);
472   s2pp = (int*)PyArray_DATA(py_s2pp_map);
473   p2s = (int*)PyArray_DATA(py_p2s_map);
474   nsym_list = (int*)PyArray_DATA(py_nsym_list);
475   n_patom = PyArray_DIMS(py_fc)[0];
476   n_satom = PyArray_DIMS(py_fc)[1];
477 
478   phpy_perm_trans_symmetrize_compact_fc(
479     fc, p2s, s2pp, nsym_list, perms, n_satom, n_patom, level);
480 
481   Py_RETURN_NONE;
482 }
483 
py_transpose_compact_fc(PyObject * self,PyObject * args)484 static PyObject * py_transpose_compact_fc(PyObject *self, PyObject *args)
485 {
486   PyArrayObject* py_fc;
487   PyArrayObject* py_permutations;
488   PyArrayObject* py_s2pp_map;
489   PyArrayObject* py_p2s_map;
490   PyArrayObject* py_nsym_list;
491   double *fc;
492   int *s2pp;
493   int *p2s;
494   int *nsym_list;
495   int *perms;
496   int n_patom, n_satom;
497 
498   if (!PyArg_ParseTuple(args, "OOOOO",
499                         &py_fc,
500                         &py_permutations,
501                         &py_s2pp_map,
502                         &py_p2s_map,
503                         &py_nsym_list)) {
504     return NULL;
505   }
506 
507   fc = (double*)PyArray_DATA(py_fc);
508   perms = (int*)PyArray_DATA(py_permutations);
509   s2pp = (int*)PyArray_DATA(py_s2pp_map);
510   p2s = (int*)PyArray_DATA(py_p2s_map);
511   nsym_list = (int*)PyArray_DATA(py_nsym_list);
512   n_patom = PyArray_DIMS(py_fc)[0];
513   n_satom = PyArray_DIMS(py_fc)[1];
514 
515   phpy_set_index_permutation_symmetry_compact_fc(fc,
516                                                  p2s,
517                                                  s2pp,
518                                                  nsym_list,
519                                                  perms,
520                                                  n_satom,
521                                                  n_patom,
522                                                  1);
523 
524   Py_RETURN_NONE;
525 }
526 
py_get_dynamical_matrix(PyObject * self,PyObject * args)527 static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
528 {
529   PyArrayObject* py_dynamical_matrix;
530   PyArrayObject* py_force_constants;
531   PyArrayObject* py_svecs;
532   PyArrayObject* py_q;
533   PyArrayObject* py_multi;
534   PyArrayObject* py_masses;
535   PyArrayObject* py_s2p_map;
536   PyArrayObject* py_p2s_map;
537 
538   double* dm;
539   double* fc;
540   double* q;
541   double (*svecs)[3];
542   double* m;
543   long (*multi)[2];
544   long* s2p_map;
545   long* p2s_map;
546   long num_patom;
547   long num_satom;
548 
549   if (!PyArg_ParseTuple(args, "OOOOOOOO",
550                         &py_dynamical_matrix,
551                         &py_force_constants,
552                         &py_q,
553                         &py_svecs,
554                         &py_multi,
555                         &py_masses,
556                         &py_s2p_map,
557                         &py_p2s_map)) {
558     return NULL;
559   }
560 
561   dm = (double*)PyArray_DATA(py_dynamical_matrix);
562   fc = (double*)PyArray_DATA(py_force_constants);
563   q = (double*)PyArray_DATA(py_q);
564   svecs = (double(*)[3])PyArray_DATA(py_svecs);
565   m = (double*)PyArray_DATA(py_masses);
566   multi = (long(*)[2])PyArray_DATA(py_multi);
567   s2p_map = (long*)PyArray_DATA(py_s2p_map);
568   p2s_map = (long*)PyArray_DATA(py_p2s_map);
569   num_patom = PyArray_DIMS(py_p2s_map)[0];
570   num_satom = PyArray_DIMS(py_s2p_map)[0];
571 
572   phpy_get_dynamical_matrix_at_q(dm,
573                                  num_patom,
574                                  num_satom,
575                                  fc,
576                                  q,
577                                  svecs,
578                                  multi,
579                                  m,
580                                  s2p_map,
581                                  p2s_map,
582                                  NULL,
583                                  1);
584 
585   Py_RETURN_NONE;
586 }
587 
588 
py_get_nac_dynamical_matrix(PyObject * self,PyObject * args)589 static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
590 {
591   PyArrayObject* py_dynamical_matrix;
592   PyArrayObject* py_force_constants;
593   PyArrayObject* py_svecs;
594   PyArrayObject* py_q_cart;
595   PyArrayObject* py_q;
596   PyArrayObject* py_multi;
597   PyArrayObject* py_masses;
598   PyArrayObject* py_s2p_map;
599   PyArrayObject* py_p2s_map;
600   PyArrayObject* py_born;
601   double factor;
602 
603   double* dm;
604   double* fc;
605   double* q_cart;
606   double* q;
607   double (*svecs)[3];
608   double* m;
609   double (*born)[3][3];
610   long (*multi)[2];
611   long* s2p_map;
612   long* p2s_map;
613   long num_patom;
614   long num_satom;
615 
616   long n;
617   double (*charge_sum)[3][3];
618 
619   if (!PyArg_ParseTuple(args, "OOOOOOOOOOd",
620                         &py_dynamical_matrix,
621                         &py_force_constants,
622                         &py_q,
623                         &py_svecs,
624                         &py_multi,
625                         &py_masses,
626                         &py_s2p_map,
627                         &py_p2s_map,
628                         &py_q_cart,
629                         &py_born,
630                         &factor))
631     return NULL;
632 
633   dm = (double*)PyArray_DATA(py_dynamical_matrix);
634   fc = (double*)PyArray_DATA(py_force_constants);
635   q_cart = (double*)PyArray_DATA(py_q_cart);
636   q = (double*)PyArray_DATA(py_q);
637   svecs = (double(*)[3])PyArray_DATA(py_svecs);
638   m = (double*)PyArray_DATA(py_masses);
639   born = (double(*)[3][3])PyArray_DATA(py_born);
640   multi = (long(*)[2])PyArray_DATA(py_multi);
641   s2p_map = (long*)PyArray_DATA(py_s2p_map);
642   p2s_map = (long*)PyArray_DATA(py_p2s_map);
643   num_patom = PyArray_DIMS(py_p2s_map)[0];
644   num_satom = PyArray_DIMS(py_s2p_map)[0];
645 
646   charge_sum = (double(*)[3][3])
647     malloc(sizeof(double[3][3]) * num_patom * num_patom);
648   n = num_satom / num_patom;
649 
650   phpy_get_charge_sum(charge_sum, num_patom, factor / n, q_cart, born);
651   phpy_get_dynamical_matrix_at_q(dm,
652                                  num_patom,
653                                  num_satom,
654                                  fc,
655                                  q,
656                                  svecs,
657                                  multi,
658                                  m,
659                                  s2p_map,
660                                  p2s_map,
661                                  charge_sum,
662                                  1);
663 
664   free(charge_sum);
665 
666   Py_RETURN_NONE;
667 }
668 
py_get_recip_dipole_dipole(PyObject * self,PyObject * args)669 static PyObject * py_get_recip_dipole_dipole(PyObject *self, PyObject *args)
670 {
671   PyArrayObject* py_dd;
672   PyArrayObject* py_dd_q0;
673   PyArrayObject* py_G_list;
674   PyArrayObject* py_q_cart;
675   PyArrayObject* py_q_direction;
676   PyArrayObject* py_born;
677   PyArrayObject* py_dielectric;
678   PyArrayObject* py_positions;
679   double factor;
680   double lambda;
681   double tolerance;
682 
683   double* dd;
684   double* dd_q0;
685   double (*G_list)[3];
686   double* q_vector;
687   double* q_direction;
688   double (*born)[3][3];
689   double (*dielectric)[3];
690   double (*pos)[3];
691   long num_patom, num_G;
692 
693   if (!PyArg_ParseTuple(args, "OOOOOOOOddd",
694                         &py_dd,
695                         &py_dd_q0,
696                         &py_G_list,
697                         &py_q_cart,
698                         &py_q_direction,
699                         &py_born,
700                         &py_dielectric,
701                         &py_positions,
702                         &factor,
703                         &lambda,
704                         &tolerance))
705     return NULL;
706 
707 
708   dd = (double*)PyArray_DATA(py_dd);
709   dd_q0 = (double*)PyArray_DATA(py_dd_q0);
710   G_list = (double(*)[3])PyArray_DATA(py_G_list);
711   if ((PyObject*)py_q_direction == Py_None) {
712     q_direction = NULL;
713   } else {
714     q_direction = (double*)PyArray_DATA(py_q_direction);
715   }
716   q_vector = (double*)PyArray_DATA(py_q_cart);
717   born = (double(*)[3][3])PyArray_DATA(py_born);
718   dielectric = (double(*)[3])PyArray_DATA(py_dielectric);
719   pos = (double(*)[3])PyArray_DATA(py_positions);
720   num_G = PyArray_DIMS(py_G_list)[0];
721   num_patom = PyArray_DIMS(py_positions)[0];
722 
723   phpy_get_recip_dipole_dipole(dd, /* [natom, 3, natom, 3, (real, imag)] */
724                                dd_q0, /* [natom, 3, 3, (real, imag)] */
725                                G_list, /* [num_kvec, 3] */
726                                num_G,
727                                num_patom,
728                                q_vector,
729                                q_direction,
730                                born,
731                                dielectric,
732                                pos, /* [natom, 3] */
733                                factor, /* 4pi/V*unit-conv */
734                                lambda, /* 4 * Lambda^2 */
735                                tolerance);
736 
737   Py_RETURN_NONE;
738 }
739 
py_get_recip_dipole_dipole_q0(PyObject * self,PyObject * args)740 static PyObject * py_get_recip_dipole_dipole_q0(PyObject *self, PyObject *args)
741 {
742   PyArrayObject* py_dd_q0;
743   PyArrayObject* py_G_list;
744   PyArrayObject* py_born;
745   PyArrayObject* py_dielectric;
746   PyArrayObject* py_positions;
747   double lambda;
748   double tolerance;
749 
750   double* dd_q0;
751   double (*G_list)[3];
752   double (*born)[3][3];
753   double (*dielectric)[3];
754   double (*pos)[3];
755   long num_patom, num_G;
756 
757   if (!PyArg_ParseTuple(args, "OOOOOdd",
758                         &py_dd_q0,
759                         &py_G_list,
760                         &py_born,
761                         &py_dielectric,
762                         &py_positions,
763                         &lambda,
764                         &tolerance))
765     return NULL;
766 
767 
768   dd_q0 = (double*)PyArray_DATA(py_dd_q0);
769   G_list = (double(*)[3])PyArray_DATA(py_G_list);
770   born = (double(*)[3][3])PyArray_DATA(py_born);
771   dielectric = (double(*)[3])PyArray_DATA(py_dielectric);
772   pos = (double(*)[3])PyArray_DATA(py_positions);
773   num_G = PyArray_DIMS(py_G_list)[0];
774   num_patom = PyArray_DIMS(py_positions)[0];
775 
776   phpy_get_recip_dipole_dipole_q0(dd_q0, /* [natom, 3, 3, (real, imag)] */
777                                   G_list, /* [num_kvec, 3] */
778                                   num_G,
779                                   num_patom,
780                                   born,
781                                   dielectric,
782                                   pos, /* [natom, 3] */
783                                   lambda, /* 4 * Lambda^2 */
784                                   tolerance);
785 
786   Py_RETURN_NONE;
787 }
788 
py_get_derivative_dynmat(PyObject * self,PyObject * args)789 static PyObject * py_get_derivative_dynmat(PyObject *self, PyObject *args)
790 {
791   PyArrayObject* py_derivative_dynmat;
792   PyArrayObject* py_force_constants;
793   PyArrayObject* py_svecs;
794   PyArrayObject* py_lattice;
795   PyArrayObject* py_q_vector;
796   PyArrayObject* py_multi;
797   PyArrayObject* py_masses;
798   PyArrayObject* py_s2p_map;
799   PyArrayObject* py_p2s_map;
800   PyArrayObject* py_born;
801   PyArrayObject* py_dielectric;
802   PyArrayObject* py_q_direction;
803   double nac_factor;
804 
805   double* ddm;
806   double* fc;
807   double* q_vector;
808   double* lat;
809   double (*svecs)[3];
810   double* masses;
811   long (*multi)[2];
812   long* s2p_map;
813   long* p2s_map;
814   long num_patom;
815   long num_satom;
816 
817   double *born;
818   double *epsilon;
819   double *q_dir;
820 
821   if (!PyArg_ParseTuple(args, "OOOOOOOOOdOOO",
822                         &py_derivative_dynmat,
823                         &py_force_constants,
824                         &py_q_vector,
825                         &py_lattice, /* column vectors */
826                         &py_svecs,
827                         &py_multi,
828                         &py_masses,
829                         &py_s2p_map,
830                         &py_p2s_map,
831                         &nac_factor,
832                         &py_born,
833                         &py_dielectric,
834                         &py_q_direction)) {
835     return NULL;
836   }
837 
838   ddm = (double*)PyArray_DATA(py_derivative_dynmat);
839   fc = (double*)PyArray_DATA(py_force_constants);
840   q_vector = (double*)PyArray_DATA(py_q_vector);
841   lat = (double*)PyArray_DATA(py_lattice);
842   svecs = (double(*)[3])PyArray_DATA(py_svecs);
843   masses = (double*)PyArray_DATA(py_masses);
844   multi = (long(*)[2])PyArray_DATA(py_multi);
845   s2p_map = (long*)PyArray_DATA(py_s2p_map);
846   p2s_map = (long*)PyArray_DATA(py_p2s_map);
847   num_patom = PyArray_DIMS(py_p2s_map)[0];
848   num_satom = PyArray_DIMS(py_s2p_map)[0];
849 
850   if ((PyObject*)py_born == Py_None) {
851     born = NULL;
852   } else {
853     born = (double*)PyArray_DATA(py_born);
854   }
855   if ((PyObject*)py_dielectric == Py_None) {
856     epsilon = NULL;
857   } else {
858     epsilon = (double*)PyArray_DATA(py_dielectric);
859   }
860   if ((PyObject*)py_q_direction == Py_None) {
861     q_dir = NULL;
862   } else {
863     q_dir = (double*)PyArray_DATA(py_q_direction);
864   }
865 
866   phpy_get_derivative_dynmat_at_q(ddm,
867                                   num_patom,
868                                   num_satom,
869                                   fc,
870                                   q_vector,
871                                   lat,
872                                   svecs,
873                                   multi,
874                                   masses,
875                                   s2p_map,
876                                   p2s_map,
877                                   nac_factor,
878                                   born,
879                                   epsilon,
880                                   q_dir);
881 
882   Py_RETURN_NONE;
883 }
884 
885 /* Thermal properties */
py_get_thermal_properties(PyObject * self,PyObject * args)886 static PyObject * py_get_thermal_properties(PyObject *self, PyObject *args)
887 {
888   PyArrayObject* py_thermal_props;
889   PyArrayObject* py_temperatures;
890   PyArrayObject* py_frequencies;
891   PyArrayObject* py_weights;
892 
893   double cutoff_frequency;
894 
895   double *temperatures;
896   double* freqs;
897   double *thermal_props;
898   long* weights;
899   long num_qpoints;
900   long num_bands;
901   long num_temp;
902 
903   if (!PyArg_ParseTuple(args, "OOOOd",
904                         &py_thermal_props,
905                         &py_temperatures,
906                         &py_frequencies,
907                         &py_weights,
908                         &cutoff_frequency)) {
909     return NULL;
910   }
911 
912   thermal_props = (double*)PyArray_DATA(py_thermal_props);
913   temperatures = (double*)PyArray_DATA(py_temperatures);
914   num_temp = (long)PyArray_DIMS(py_temperatures)[0];
915   freqs = (double*)PyArray_DATA(py_frequencies);
916   num_qpoints = (long)PyArray_DIMS(py_frequencies)[0];
917   weights = (long*)PyArray_DATA(py_weights);
918   num_bands = (long)PyArray_DIMS(py_frequencies)[1];
919 
920   phpy_get_thermal_properties(thermal_props,
921                               temperatures,
922                               freqs,
923                               weights,
924                               num_temp,
925                               num_qpoints,
926                               num_bands,
927                               cutoff_frequency);
928 
929   Py_RETURN_NONE;
930 }
931 
py_distribute_fc2(PyObject * self,PyObject * args)932 static PyObject * py_distribute_fc2(PyObject *self, PyObject *args)
933 {
934   PyArrayObject* py_force_constants;
935   PyArrayObject* py_permutations;
936   PyArrayObject* py_map_atoms;
937   PyArrayObject* py_map_syms;
938   PyArrayObject* py_atom_list;
939   PyArrayObject* py_rotations_cart;
940 
941   double (*r_carts)[3][3];
942   double (*fc2)[3][3];
943   int *permutations;
944   int *map_atoms;
945   int *map_syms;
946   int *atom_list;
947   npy_intp num_pos, num_rot, len_atom_list;
948 
949   if (!PyArg_ParseTuple(args, "OOOOOO",
950                         &py_force_constants,
951                         &py_atom_list,
952                         &py_rotations_cart,
953                         &py_permutations,
954                         &py_map_atoms,
955                         &py_map_syms)) {
956     return NULL;
957   }
958 
959   fc2 = (double(*)[3][3])PyArray_DATA(py_force_constants);
960   atom_list = (int*)PyArray_DATA(py_atom_list);
961   len_atom_list = PyArray_DIMS(py_atom_list)[0];
962   permutations = (int*)PyArray_DATA(py_permutations);
963   map_atoms = (int*)PyArray_DATA(py_map_atoms);
964   map_syms = (int*)PyArray_DATA(py_map_syms);
965   r_carts = (double(*)[3][3])PyArray_DATA(py_rotations_cart);
966   num_rot = PyArray_DIMS(py_permutations)[0];
967   num_pos = PyArray_DIMS(py_permutations)[1];
968 
969   if (PyArray_NDIM(py_map_atoms) != 1 || PyArray_DIMS(py_map_atoms)[0] != num_pos)
970   {
971     PyErr_SetString(PyExc_ValueError, "wrong shape for map_atoms");
972     return NULL;
973   }
974 
975   if (PyArray_NDIM(py_map_syms) != 1 || PyArray_DIMS(py_map_syms)[0] != num_pos)
976   {
977     PyErr_SetString(PyExc_ValueError, "wrong shape for map_syms");
978     return NULL;
979   }
980 
981   if (PyArray_DIMS(py_rotations_cart)[0] != num_rot)
982   {
983     PyErr_SetString(PyExc_ValueError, "permutations and rotations are different length");
984     return NULL;
985   }
986 
987   phpy_distribute_fc2(fc2,
988                       atom_list,
989                       len_atom_list,
990                       r_carts,
991                       permutations,
992                       map_atoms,
993                       map_syms,
994                       num_rot,
995                       num_pos);
996   Py_RETURN_NONE;
997 }
998 
999 
1000 static PyObject *
py_thm_relative_grid_address(PyObject * self,PyObject * args)1001 py_thm_relative_grid_address(PyObject *self, PyObject *args)
1002 {
1003   PyArrayObject* py_relative_grid_address;
1004   PyArrayObject* py_reciprocal_lattice_py;
1005 
1006   long (*relative_grid_address)[4][3];
1007   double (*reciprocal_lattice)[3];
1008 
1009   if (!PyArg_ParseTuple(args, "OO",
1010                         &py_relative_grid_address,
1011                         &py_reciprocal_lattice_py)) {
1012     return NULL;
1013   }
1014 
1015   relative_grid_address = (long(*)[4][3])PyArray_DATA(py_relative_grid_address);
1016   reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice_py);
1017 
1018   phpy_get_relative_grid_address(relative_grid_address, reciprocal_lattice);
1019 
1020   Py_RETURN_NONE;
1021 }
1022 
1023 static PyObject *
py_thm_all_relative_grid_address(PyObject * self,PyObject * args)1024 py_thm_all_relative_grid_address(PyObject *self, PyObject *args)
1025 {
1026   PyArrayObject* py_relative_grid_address;
1027 
1028   long (*relative_grid_address)[24][4][3];
1029 
1030   if (!PyArg_ParseTuple(args, "O",
1031                         &py_relative_grid_address)) {
1032     return NULL;
1033   }
1034 
1035   relative_grid_address =
1036     (long(*)[24][4][3])PyArray_DATA(py_relative_grid_address);
1037 
1038   phpy_get_all_relative_grid_address(relative_grid_address);
1039 
1040   Py_RETURN_NONE;
1041 }
1042 
1043 static PyObject *
py_thm_integration_weight(PyObject * self,PyObject * args)1044 py_thm_integration_weight(PyObject *self, PyObject *args)
1045 {
1046   double omega;
1047   PyArrayObject* py_tetrahedra_omegas;
1048   char* function;
1049 
1050   double (*tetrahedra_omegas)[4];
1051   double iw;
1052 
1053   if (!PyArg_ParseTuple(args, "dOs",
1054                         &omega,
1055                         &py_tetrahedra_omegas,
1056                         &function)) {
1057     return NULL;
1058   }
1059 
1060   tetrahedra_omegas = (double(*)[4])PyArray_DATA(py_tetrahedra_omegas);
1061 
1062   iw = phpy_get_integration_weight(omega,
1063                                    tetrahedra_omegas,
1064                                    function[0]);
1065 
1066   return PyFloat_FromDouble(iw);
1067 }
1068 
1069 static PyObject *
py_thm_integration_weight_at_omegas(PyObject * self,PyObject * args)1070 py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args)
1071 {
1072   PyArrayObject* py_integration_weights;
1073   PyArrayObject* py_omegas;
1074   PyArrayObject* py_tetrahedra_omegas;
1075   char* function;
1076 
1077   double *omegas;
1078   double *iw;
1079   long num_omegas;
1080   double (*tetrahedra_omegas)[4];
1081 
1082   long i;
1083 
1084   if (!PyArg_ParseTuple(args, "OOOs",
1085                         &py_integration_weights,
1086                         &py_omegas,
1087                         &py_tetrahedra_omegas,
1088                         &function)) {
1089     return NULL;
1090   }
1091 
1092   omegas = (double*)PyArray_DATA(py_omegas);
1093   iw = (double*)PyArray_DATA(py_integration_weights);
1094   num_omegas = (long)PyArray_DIMS(py_omegas)[0];
1095   tetrahedra_omegas = (double(*)[4])PyArray_DATA(py_tetrahedra_omegas);
1096 
1097 #pragma omp parallel for
1098   for (i = 0; i < num_omegas; i++) {
1099     iw[i] = phpy_get_integration_weight(omegas[i],
1100                                         tetrahedra_omegas,
1101                                         function[0]);
1102   }
1103 
1104   Py_RETURN_NONE;
1105 }
1106 
py_get_tetrahedra_frequenies(PyObject * self,PyObject * args)1107 static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args)
1108 {
1109   PyArrayObject* py_freq_tetras;
1110   PyArrayObject* py_grid_points;
1111   PyArrayObject* py_mesh;
1112   PyArrayObject* py_grid_address;
1113   PyArrayObject* py_gp_ir_index;
1114   PyArrayObject* py_relative_grid_address;
1115   PyArrayObject* py_frequencies;
1116 
1117   double* freq_tetras;
1118   long* grid_points;
1119   long* mesh;
1120   long (*grid_address)[3];
1121   long* gp_ir_index;
1122   long (*relative_grid_address)[3];
1123   double* frequencies;
1124 
1125   long num_gp_in, num_band;
1126 
1127   if (!PyArg_ParseTuple(args, "OOOOOOO",
1128                         &py_freq_tetras,
1129                         &py_grid_points,
1130                         &py_mesh,
1131                         &py_grid_address,
1132                         &py_gp_ir_index,
1133                         &py_relative_grid_address,
1134                         &py_frequencies)) {
1135     return NULL;
1136   }
1137 
1138   freq_tetras = (double*)PyArray_DATA(py_freq_tetras);
1139   grid_points = (long*)PyArray_DATA(py_grid_points);
1140   num_gp_in = PyArray_DIMS(py_grid_points)[0];
1141   mesh = (long*)PyArray_DATA(py_mesh);
1142   grid_address = (long(*)[3])PyArray_DATA(py_grid_address);
1143   gp_ir_index = (long*)PyArray_DATA(py_gp_ir_index);
1144   relative_grid_address = (long(*)[3])PyArray_DATA(py_relative_grid_address);
1145   frequencies = (double*)PyArray_DATA(py_frequencies);
1146   num_band = PyArray_DIMS(py_frequencies)[1];
1147 
1148   phpy_get_tetrahedra_frequenies(freq_tetras,
1149                                  mesh,
1150                                  grid_points,
1151                                  grid_address,
1152                                  relative_grid_address,
1153                                  gp_ir_index,
1154                                  frequencies,
1155                                  num_band,
1156                                  num_gp_in);
1157 
1158   Py_RETURN_NONE;
1159 }
1160 
py_tetrahedron_method_dos(PyObject * self,PyObject * args)1161 static PyObject * py_tetrahedron_method_dos(PyObject *self, PyObject *args)
1162 {
1163   PyArrayObject* py_dos;
1164   PyArrayObject* py_mesh;
1165   PyArrayObject* py_freq_points;
1166   PyArrayObject* py_frequencies;
1167   PyArrayObject* py_coef;
1168   PyArrayObject* py_grid_address;
1169   PyArrayObject* py_grid_mapping_table;
1170   PyArrayObject* py_relative_grid_address;
1171 
1172   double *dos;
1173   long* mesh;
1174   double* freq_points;
1175   double* frequencies;
1176   double* coef;
1177   long (*grid_address)[3];
1178   long num_gp, num_ir_gp, num_band, num_freq_points, num_coef;
1179   long *grid_mapping_table;
1180   long (*relative_grid_address)[4][3];
1181 
1182   if (!PyArg_ParseTuple(args, "OOOOOOOO",
1183                         &py_dos,
1184                         &py_mesh,
1185                         &py_freq_points,
1186                         &py_frequencies,
1187                         &py_coef,
1188                         &py_grid_address,
1189                         &py_grid_mapping_table,
1190                         &py_relative_grid_address)) {
1191     return NULL;
1192   }
1193 
1194   /* dos[num_ir_gp][num_band][num_freq_points][num_coef] */
1195   dos = (double*)PyArray_DATA(py_dos);
1196   mesh = (long*)PyArray_DATA(py_mesh);
1197   freq_points = (double*)PyArray_DATA(py_freq_points);
1198   num_freq_points = (long)PyArray_DIMS(py_freq_points)[0];
1199   frequencies = (double*)PyArray_DATA(py_frequencies);
1200   num_ir_gp = (long)PyArray_DIMS(py_frequencies)[0];
1201   num_band = (long)PyArray_DIMS(py_frequencies)[1];
1202   coef = (double*)PyArray_DATA(py_coef);
1203   num_coef = (long)PyArray_DIMS(py_coef)[1];
1204   grid_address = (long(*)[3])PyArray_DATA(py_grid_address);
1205   num_gp = (long)PyArray_DIMS(py_grid_address)[0];
1206   grid_mapping_table = (long*)PyArray_DATA(py_grid_mapping_table);
1207   relative_grid_address = (long(*)[4][3])PyArray_DATA(py_relative_grid_address);
1208 
1209   phpy_tetrahedron_method_dos(dos,
1210                               mesh,
1211                               grid_address,
1212                               relative_grid_address,
1213                               grid_mapping_table,
1214                               freq_points,
1215                               frequencies,
1216                               coef,
1217                               num_freq_points,
1218                               num_ir_gp,
1219                               num_band,
1220                               num_coef,
1221                               num_gp);
1222 
1223   Py_RETURN_NONE;
1224 }
1225