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