1 /* Copyright (C) 2010-2011 CAMd
2 * Please see the accompanying LICENSE file for further information. */
3 #include "extensions.h"
4
5 //
6 // Apply symmetry operation op_cc to a and add result to b:
7 //
8 // =T_ _
9 // b(U g) += a(g),
10 //
11 // where:
12 //
13 // = _T
14 // U = op_cc[c1, c2] and g = (g0, g1, g2).
15 // c1,c2
16 //
symmetrize(PyObject * self,PyObject * args)17 PyObject* symmetrize(PyObject *self, PyObject *args)
18 {
19 PyArrayObject* a_g_obj;
20 PyArrayObject* b_g_obj;
21 PyArrayObject* op_cc_obj;
22 PyArrayObject* offset_c_obj;
23
24 if (!PyArg_ParseTuple(args, "OOOO",
25 &a_g_obj, &b_g_obj, &op_cc_obj, &offset_c_obj))
26 return NULL;
27
28 const long* C = (const long*)PyArray_DATA(op_cc_obj);
29 const long* o_c = (const long*)PyArray_DATA(offset_c_obj);
30 int ng0 = PyArray_DIMS(a_g_obj)[0];
31 int ng1 = PyArray_DIMS(a_g_obj)[1];
32 int ng2 = PyArray_DIMS(a_g_obj)[2];
33 int Ng0 = ng0 + o_c[0];
34 int Ng1 = ng1 + o_c[1];
35 int Ng2 = ng2 + o_c[2];
36
37 const double* a_g = (const double*)PyArray_DATA(a_g_obj);
38 double* b_g = (double*)PyArray_DATA(b_g_obj);
39 for (int g0 = o_c[0]; g0 < Ng0; g0++)
40 for (int g1 = o_c[1]; g1 < Ng1; g1++)
41 for (int g2 = o_c[2]; g2 < Ng2; g2++) {
42 int p0 = ((C[0] * g0 + C[3] * g1 + C[6] * g2) %
43 Ng0 + Ng0) % Ng0;
44 int p1 = ((C[1] * g0 + C[4] * g1 + C[7] * g2) %
45 Ng1 + Ng1) % Ng1;
46 int p2 = ((C[2] * g0 + C[5] * g1 + C[8] * g2) %
47 Ng2 + Ng2) % Ng2;
48 b_g[((p0 - o_c[0]) * ng1 +
49 (p1 - o_c[1])) * ng2 +
50 p2 - o_c[2]] += *a_g++;
51 }
52
53 Py_RETURN_NONE;
54 }
55
symmetrize_ft(PyObject * self,PyObject * args)56 PyObject* symmetrize_ft(PyObject *self, PyObject *args)
57 {
58 PyArrayObject* a_g_obj;
59 PyArrayObject* b_g_obj;
60 PyArrayObject* op_cc_obj;
61 PyArrayObject* t_c_obj;
62 PyArrayObject* offset_c_obj;
63
64 if (!PyArg_ParseTuple(args, "OOOOO",
65 &a_g_obj, &b_g_obj, &op_cc_obj, &t_c_obj,
66 &offset_c_obj))
67 return NULL;
68
69 const long* t_c = (const long*)PyArray_DATA(t_c_obj);
70 const long* C = (const long*)PyArray_DATA(op_cc_obj);
71 const long* o_c = (const long*)PyArray_DATA(offset_c_obj);
72
73 int ng0 = PyArray_DIMS(a_g_obj)[0];
74 int ng1 = PyArray_DIMS(a_g_obj)[1];
75 int ng2 = PyArray_DIMS(a_g_obj)[2];
76 int Ng0 = ng0 + o_c[0];
77 int Ng1 = ng1 + o_c[1];
78 int Ng2 = ng2 + o_c[2];
79
80 const double* a_g = (const double*)PyArray_DATA(a_g_obj);
81 double* b_g = (double*)PyArray_DATA(b_g_obj);
82 for (int g0 = 0; g0 < ng0; g0++)
83 for (int g1 = 0; g1 < ng1; g1++)
84 for (int g2 = 0; g2 < ng2; g2++) {
85 int p0 = ((C[0] * g0 + C[3] * g1 + C[6] * g2 - t_c[0]) %
86 Ng0 + Ng0) % Ng0;
87 int p1 = ((C[1] * g0 + C[4] * g1 + C[7] * g2 - t_c[1]) %
88 Ng1 + Ng1) % Ng1;
89 int p2 = ((C[2] * g0 + C[5] * g1 + C[8] * g2 - t_c[2]) %
90 Ng2 + Ng2) % Ng2;
91 b_g[((p0 - o_c[0]) * ng1 +
92 (p1 - o_c[1])) * ng2 +
93 p2 - o_c[2]] += *a_g++;
94 }
95
96 Py_RETURN_NONE;
97 }
98
symmetrize_wavefunction(PyObject * self,PyObject * args)99 PyObject* symmetrize_wavefunction(PyObject *self, PyObject *args)
100 {
101 PyArrayObject* a_g_obj;
102 PyArrayObject* b_g_obj;
103 PyArrayObject* op_cc_obj;
104 PyArrayObject* kpt0_obj;
105 PyArrayObject* kpt1_obj;
106
107 if (!PyArg_ParseTuple(args, "OOOOO", &a_g_obj, &b_g_obj, &op_cc_obj, &kpt0_obj, &kpt1_obj))
108 return NULL;
109
110 const long* C = (const long*)PyArray_DATA(op_cc_obj);
111 const double* kpt0 = (const double*) PyArray_DATA(kpt0_obj);
112 const double* kpt1 = (const double*) PyArray_DATA(kpt1_obj);
113 int ng0 = PyArray_DIMS(a_g_obj)[0];
114 int ng1 = PyArray_DIMS(a_g_obj)[1];
115 int ng2 = PyArray_DIMS(a_g_obj)[2];
116
117 const double complex* a_g = (const double complex*)PyArray_DATA(a_g_obj);
118 double complex* b_g = (double complex*)PyArray_DATA(b_g_obj);
119
120 for (int g0 = 0; g0 < ng0; g0++)
121 for (int g1 = 0; g1 < ng1; g1++)
122 for (int g2 = 0; g2 < ng2; g2++) {
123 int p0 = ((C[0] * g0 + C[3] * g1 + C[6] * g2) % ng0 + ng0) % ng0;
124 int p1 = ((C[1] * g0 + C[4] * g1 + C[7] * g2) % ng1 + ng1) % ng1;
125 int p2 = ((C[2] * g0 + C[5] * g1 + C[8] * g2) % ng2 + ng2) % ng2;
126
127 double complex phase = cexp(I * 2. * M_PI *
128 (kpt1[0]/ng0*p0 +
129 kpt1[1]/ng1*p1 +
130 kpt1[2]/ng2*p2 -
131 kpt0[0]/ng0*g0 -
132 kpt0[1]/ng1*g1 -
133 kpt0[2]/ng2*g2));
134 b_g[(p0 * ng1 + p1) * ng2 + p2] += (*a_g * phase);
135 a_g++;
136 }
137
138 Py_RETURN_NONE;
139 }
140
symmetrize_return_index(PyObject * self,PyObject * args)141 PyObject* symmetrize_return_index(PyObject *self, PyObject *args)
142 {
143 PyArrayObject* a_g_obj;
144 PyArrayObject* b_g_obj;
145 PyArrayObject* op_cc_obj;
146 PyArrayObject* kpt0_obj;
147 PyArrayObject* kpt1_obj;
148
149 if (!PyArg_ParseTuple(args, "OOOOO", &a_g_obj, &b_g_obj, &op_cc_obj, &kpt0_obj, &kpt1_obj))
150 return NULL;
151
152 const long* C = (const long*)PyArray_DATA(op_cc_obj);
153 const double* kpt0 = (const double*) PyArray_DATA(kpt0_obj);
154 const double* kpt1 = (const double*) PyArray_DATA(kpt1_obj);
155
156 int ng0 = PyArray_DIMS(a_g_obj)[0];
157 int ng1 = PyArray_DIMS(a_g_obj)[1];
158 int ng2 = PyArray_DIMS(a_g_obj)[2];
159
160 unsigned long* a_g = (unsigned long*)PyArray_DATA(a_g_obj);
161 double complex* b_g = (double complex*)PyArray_DATA(b_g_obj);
162
163 for (int g0 = 0; g0 < ng0; g0++)
164 for (int g1 = 0; g1 < ng1; g1++)
165 for (int g2 = 0; g2 < ng2; g2++) {
166 int p0 = ((C[0] * g0 + C[3] * g1 + C[6] * g2) % ng0 + ng0) % ng0;
167 int p1 = ((C[1] * g0 + C[4] * g1 + C[7] * g2) % ng1 + ng1) % ng1;
168 int p2 = ((C[2] * g0 + C[5] * g1 + C[8] * g2) % ng2 + ng2) % ng2;
169
170 double complex phase = cexp(I * 2. * M_PI *
171 (kpt1[0]/ng0*p0 +
172 kpt1[1]/ng1*p1 +
173 kpt1[2]/ng2*p2 -
174 kpt0[0]/ng0*g0 -
175 kpt0[1]/ng1*g1 -
176 kpt0[2]/ng2*g2));
177 *a_g++ = (p0 * ng1 + p1) * ng2 + p2;
178 *b_g++ = phase;
179 }
180
181 Py_RETURN_NONE;
182 }
183
symmetrize_with_index(PyObject * self,PyObject * args)184 PyObject* symmetrize_with_index(PyObject *self, PyObject *args)
185 {
186 PyArrayObject* a_g_obj;
187 PyArrayObject* b_g_obj;
188 PyArrayObject* index_g_obj;
189 PyArrayObject* phase_g_obj;
190
191 if (!PyArg_ParseTuple(args, "OOOO", &a_g_obj, &b_g_obj, &index_g_obj, &phase_g_obj))
192 return NULL;
193
194 int ng0 = PyArray_DIMS(a_g_obj)[0];
195 int ng1 = PyArray_DIMS(a_g_obj)[1];
196 int ng2 = PyArray_DIMS(a_g_obj)[2];
197
198 const unsigned long* index_g = (const unsigned long*)PyArray_DATA(index_g_obj);
199 const double complex* phase_g = (const double complex*)PyArray_DATA(phase_g_obj);
200 const double complex* a_g = (const double complex*)PyArray_DATA(a_g_obj);
201 double complex* b_g = (double complex*)PyArray_DATA(b_g_obj);
202
203
204 for (int g0 = 0; g0 < ng0; g0++)
205 for (int g1 = 0; g1 < ng1; g1++)
206 for (int g2 = 0; g2 < ng2; g2++) {
207 b_g[*index_g] += (*a_g * *phase_g);
208 a_g++;
209 phase_g++;
210 index_g++;
211 }
212
213 Py_RETURN_NONE;
214 }
215
map_k_points(PyObject * self,PyObject * args)216 PyObject* map_k_points(PyObject *self, PyObject *args)
217 {
218 PyArrayObject* bzk_kc_obj;
219 PyArrayObject* U_scc_obj;
220 double tol;
221 PyArrayObject* bz2bz_ks_obj;
222 int ka, kb;
223
224 if (!PyArg_ParseTuple(args, "OOdOii", &bzk_kc_obj, &U_scc_obj,
225 &tol, &bz2bz_ks_obj, &ka, &kb))
226 return NULL;
227
228 const long* U_scc = (const long*)PyArray_DATA(U_scc_obj);
229 const double* bzk_kc = (const double*)PyArray_DATA(bzk_kc_obj);
230 long* bz2bz_ks = (long*)PyArray_DATA(bz2bz_ks_obj);
231
232 int nbzkpts = PyArray_DIMS(bzk_kc_obj)[0];
233 int nsym = PyArray_DIMS(U_scc_obj)[0];
234
235 for (int k1 = ka; k1 < kb; k1++) {
236 const double* q = bzk_kc + k1 * 3;
237 for (int s = 0; s < nsym; s++) {
238 const long* U = U_scc + s * 9;
239 double q0 = U[0] * q[0] + U[1] * q[1] + U[2] * q[2];
240 double q1 = U[3] * q[0] + U[4] * q[1] + U[5] * q[2];
241 double q2 = U[6] * q[0] + U[7] * q[1] + U[8] * q[2];
242 for (int k2 = 0; k2 < nbzkpts; k2++) {
243 double p0 = q0 - bzk_kc[k2 * 3];
244 if (fabs(p0 - round(p0)) > tol)
245 continue;
246 double p1 = q1 - bzk_kc[k2 * 3 + 1];
247 if (fabs(p1 - round(p1)) > tol)
248 continue;
249 double p2 = q2 - bzk_kc[k2 * 3 + 2];
250 if (fabs(p2 - round(p2)) > tol)
251 continue;
252 bz2bz_ks[k1 * nsym + s] = k2;
253 break;
254 }
255 }
256 }
257 Py_RETURN_NONE;
258 }
259