1 /* Copyright (C) 2010 CAMd
2 * Please see the accompanying LICENSE file for further information. */
3
4 #include "extensions.h"
5 #include "spline.h"
6 #include "lfc.h"
7 #include "bmgs/spherical_harmonics.h"
8
second_derivative(LFCObject * lfc,PyObject * args)9 PyObject* second_derivative(LFCObject *lfc, PyObject *args)
10 {
11 PyArrayObject* a_G_obj;
12 PyArrayObject* c_Mvv_obj;
13 PyArrayObject* h_cv_obj;
14 PyArrayObject* n_c_obj;
15 PyObject* spline_M_obj;
16 PyArrayObject* beg_c_obj;
17 PyArrayObject* pos_Wc_obj;
18 int q;
19
20 if (!PyArg_ParseTuple(args, "OOOOOOOi", &a_G_obj, &c_Mvv_obj,
21 &h_cv_obj, &n_c_obj,
22 &spline_M_obj, &beg_c_obj,
23 &pos_Wc_obj, &q))
24 return NULL;
25
26 // Copied from derivative member function
27 int nd = PyArray_NDIM(a_G_obj);
28 npy_intp* dims = PyArray_DIMS(a_G_obj);
29 int nx = PyArray_MultiplyList(dims, nd - 3);
30 int nG = PyArray_MultiplyList(dims + nd - 3, 3);
31 int nM = PyArray_DIM(c_Mvv_obj, PyArray_NDIM(c_Mvv_obj) - 2);
32
33 // These were already present
34 const double* h_cv = (const double*)PyArray_DATA(h_cv_obj);
35 const long* n_c = (const long*)PyArray_DATA(n_c_obj);
36 const double (*pos_Wc)[3] = (const double (*)[3])PyArray_DATA(pos_Wc_obj);
37
38 long* beg_c = LONGP(beg_c_obj);
39 ///////////////////////////////////////////////
40
41 const double Y00dv = lfc->dv / sqrt(4.0 * M_PI);
42
43 if (!lfc->bloch_boundary_conditions) {
44 const double* a_G = (const double*)PyArray_DATA(a_G_obj);
45 double* c_Mvv = (double*)PyArray_DATA(c_Mvv_obj);
46 // Loop over number of x-dimension in a_xG (not relevant yet)
47 for (int x = 0; x < nx; x++) {
48 // JJs old stuff
49 GRID_LOOP_START(lfc, -1, 0) {
50 // In one grid loop iteration, only i2 changes.
51 int i2 = Ga % n_c[2] + beg_c[2];
52 int i1 = (Ga / n_c[2]) % n_c[1] + beg_c[1];
53 int i0 = Ga / (n_c[2] * n_c[1]) + beg_c[0];
54 double xG = h_cv[0] * i0 + h_cv[3] * i1 + h_cv[6] * i2;
55 double yG = h_cv[1] * i0 + h_cv[4] * i1 + h_cv[7] * i2;
56 double zG = h_cv[2] * i0 + h_cv[5] * i1 + h_cv[8] * i2;
57 for (int G = Ga; G < Gb; G++) {
58 for (int i = 0; i < ni; i++) {
59 LFVolume* vol = volume_i[i];
60 int M = vol->M;
61 double* c_mvv = c_Mvv + 9 * M;
62 const bmgsspline* spline = (const bmgsspline*) \
63 &((const SplineObject*)PyList_GetItem(spline_M_obj, M))->spline;
64
65 double x = xG - pos_Wc[vol->W][0];
66 double y = yG - pos_Wc[vol->W][1];
67 double z = zG - pos_Wc[vol->W][2];
68 double r2 = x * x + y * y + z * z;
69 double r = sqrt(r2);
70 int bin = r / spline->dr;
71 assert(bin <= spline->nbins);
72 double* s = spline->data + 4 * bin;
73 double u = r - bin * spline->dr;
74 double dfdror;
75 if (bin == 0)
76 dfdror = 2.0 * s[2] + 3.0 * s[3] * r;
77 else
78 dfdror = (s[1] + u * (2.0 * s[2] + u * 3.0 * s[3])) / r;
79 double a = a_G[G] * Y00dv;
80 dfdror *= a;
81 c_mvv[0] += dfdror;
82 c_mvv[4] += dfdror;
83 c_mvv[8] += dfdror;
84 if (r > 1e-15) {
85 double b = ((2.0 * s[2] + 6.0 * s[3] * u) * a - dfdror) / r2;
86 c_mvv[0] += b * x * x;
87 c_mvv[1] += b * x * y;
88 c_mvv[2] += b * x * z;
89 c_mvv[3] += b * y * x;
90 c_mvv[4] += b * y * y;
91 c_mvv[5] += b * y * z;
92 c_mvv[6] += b * z * x;
93 c_mvv[7] += b * z * y;
94 c_mvv[8] += b * z * z;
95 }
96 }
97 xG += h_cv[6];
98 yG += h_cv[7];
99 zG += h_cv[8];
100 }
101 }
102 GRID_LOOP_STOP(lfc, -1, 0);
103 c_Mvv += 9 * nM;
104 a_G += nG;
105 }
106 }
107 else {
108 const complex double* a_G = (const complex double*)PyArray_DATA(a_G_obj);
109 complex double* c_Mvv = (complex double*)PyArray_DATA(c_Mvv_obj);
110
111 for (int x = 0; x < nx; x++) {
112 GRID_LOOP_START(lfc, q, 0) {
113 // In one grid loop iteration, only i2 changes.
114 int i2 = Ga % n_c[2] + beg_c[2];
115 int i1 = (Ga / n_c[2]) % n_c[1] + beg_c[1];
116 int i0 = Ga / (n_c[2] * n_c[1]) + beg_c[0];
117 double xG = h_cv[0] * i0 + h_cv[3] * i1 + h_cv[6] * i2;
118 double yG = h_cv[1] * i0 + h_cv[4] * i1 + h_cv[7] * i2;
119 double zG = h_cv[2] * i0 + h_cv[5] * i1 + h_cv[8] * i2;
120 for (int G = Ga; G < Gb; G++) {
121 for (int i = 0; i < ni; i++) {
122 LFVolume* vol = volume_i[i];
123 int M = vol->M;
124 complex double* c_mvv = c_Mvv + 9 * M;
125 const bmgsspline* spline = (const bmgsspline*) \
126 &((const SplineObject*)PyList_GetItem(spline_M_obj, M))->spline;
127
128 double x = xG - pos_Wc[vol->W][0];
129 double y = yG - pos_Wc[vol->W][1];
130 double z = zG - pos_Wc[vol->W][2];
131 double r2 = x * x + y * y + z * z;
132 double r = sqrt(r2);
133 double dfdror;
134
135 // use bmgs_get_value_and_derivative instead ??!!
136 int bin = r / spline->dr;
137 assert(bin <= spline->nbins);
138 double u = r - bin * spline->dr;
139 double* s = spline->data + 4 * bin;
140
141 if (bin == 0)
142 dfdror = 2.0 * s[2] + 3.0 * s[3] * r;
143 else
144 dfdror = (s[1] + u * (2.0 * s[2] + u * 3.0 * s[3])) / r;
145 // phase added here
146 complex double a = a_G[G] * phase_i[i] * Y00dv;
147 // dfdror *= a;
148 c_mvv[0] += a * dfdror;
149 c_mvv[4] += a * dfdror;
150 c_mvv[8] += a * dfdror;
151 if (r > 1e-15) {
152 double b = (2.0 * s[2] + 6.0 * s[3] * u - dfdror) / r2;
153 c_mvv[0] += a * b * x * x;
154 c_mvv[1] += a * b * x * y;
155 c_mvv[2] += a * b * x * z;
156 c_mvv[3] += a * b * y * x;
157 c_mvv[4] += a * b * y * y;
158 c_mvv[5] += a * b * y * z;
159 c_mvv[6] += a * b * z * x;
160 c_mvv[7] += a * b * z * y;
161 c_mvv[8] += a * b * z * z;
162 }
163 }
164 xG += h_cv[6];
165 yG += h_cv[7];
166 zG += h_cv[8];
167 }
168 }
169 GRID_LOOP_STOP(lfc, q, 0);
170 c_Mvv += 9 * nM;
171 a_G += nG;
172 }
173 }
174 Py_RETURN_NONE;
175 }
176
add_derivative(LFCObject * lfc,PyObject * args)177 PyObject* add_derivative(LFCObject *lfc, PyObject *args)
178 {
179 // Coefficients for the lfc's
180 PyArrayObject* c_xM_obj;
181 // Array
182 PyArrayObject* a_xG_obj;
183 PyArrayObject* h_cv_obj;
184 PyArrayObject* n_c_obj;
185 PyObject* spline_M_obj;
186 PyArrayObject* beg_c_obj;
187 PyArrayObject* pos_Wc_obj;
188 // Atom index
189 int a;
190 // Cartesian coordinate
191 int v;
192 // k-point index
193 int q;
194
195 if (!PyArg_ParseTuple(args, "OOOOOOOiii", &c_xM_obj, &a_xG_obj,
196 &h_cv_obj, &n_c_obj, &spline_M_obj, &beg_c_obj,
197 &pos_Wc_obj, &a, &v, &q))
198 return NULL;
199
200 // Number of dimensions
201 int nd = PyArray_NDIM(a_xG_obj);
202 // Array with lengths of array dimensions
203 npy_intp* dims = PyArray_DIMS(a_xG_obj);
204 // Number of extra dimensions
205 int nx = PyArray_MultiplyList(dims, nd - 3);
206 // Number of grid points
207 int nG = PyArray_MultiplyList(dims + nd - 3, 3);
208 // Number of lfc's
209 int nM = PyArray_DIM(c_xM_obj, PyArray_NDIM(c_xM_obj) - 1);
210
211 const double* h_cv = (const double*)PyArray_DATA(h_cv_obj);
212 const long* n_c = (const long*)PyArray_DATA(n_c_obj);
213 const double (*pos_Wc)[3] = (const double (*)[3])PyArray_DATA(pos_Wc_obj);
214
215 long* beg_c = LONGP(beg_c_obj);
216
217 if (!lfc->bloch_boundary_conditions) {
218
219 const double* c_M = (const double*)PyArray_DATA(c_xM_obj);
220 double* a_G = (double*)PyArray_DATA(a_xG_obj);
221 for (int x = 0; x < nx; x++) {
222 GRID_LOOP_START(lfc, -1, 0) {
223
224 // In one grid loop iteration, only i2 changes.
225 int i2 = Ga % n_c[2] + beg_c[2];
226 int i1 = (Ga / n_c[2]) % n_c[1] + beg_c[1];
227 int i0 = Ga / (n_c[2] * n_c[1]) + beg_c[0];
228 // Grid point position
229 double xG = h_cv[0] * i0 + h_cv[3] * i1 + h_cv[6] * i2;
230 double yG = h_cv[1] * i0 + h_cv[4] * i1 + h_cv[7] * i2;
231 double zG = h_cv[2] * i0 + h_cv[5] * i1 + h_cv[8] * i2;
232 // Loop over grid points in current stride
233 for (int G = Ga; G < Gb; G++) {
234 // Loop over volumes at current grid point
235 for (int i = 0; i < ni; i++) {
236
237 LFVolume* vol = volume_i[i];
238 int M = vol->M;
239 // Check that the volume belongs to the atom in consideration later
240 int W = vol->W;
241 int nm = vol->nm;
242 int l = (nm - 1) / 2;
243
244 const bmgsspline* spline = (const bmgsspline*) \
245 &((const SplineObject*)PyList_GetItem(spline_M_obj, M))->spline;
246
247 double x = xG - pos_Wc[W][0];
248 double y = yG - pos_Wc[W][1];
249 double z = zG - pos_Wc[W][2];
250 double R_c[] = {x, y, z};
251 double r2 = x * x + y * y + z * z;
252 double r = sqrt(r2);
253 double f;
254 double dfdr;
255
256 bmgs_get_value_and_derivative(spline, r, &f, &dfdr);
257
258 // First contribution: f * d(r^l * Y)/dv
259 double fdrlYdx_m[nm];
260 if (v == 0)
261 spherical_harmonics_derivative_x(l, f, x, y, z, r2, fdrlYdx_m);
262 else if (v == 1)
263 spherical_harmonics_derivative_y(l, f, x, y, z, r2, fdrlYdx_m);
264 else
265 spherical_harmonics_derivative_z(l, f, x, y, z, r2, fdrlYdx_m);
266
267 for (int m = 0; m < nm; m++)
268 a_G[G] += fdrlYdx_m[m] * c_M[M + m];
269
270 // Second contribution: r^(l-1) * Y * df/dr * R_v
271 if (r > 1e-15) {
272 double rlm1Ydfdr_m[nm]; // r^(l-1) * Y * df/dr
273 double rm1dfdr = 1. / r * dfdr;
274 spherical_harmonics(l, rm1dfdr, x, y, z, r2, rlm1Ydfdr_m);
275 for (int m = 0; m < nm; m++)
276 a_G[G] += rlm1Ydfdr_m[m] * R_c[v] * c_M[M + m];
277 }
278 }
279 // Update coordinates of current grid point
280 xG += h_cv[6];
281 yG += h_cv[7];
282 zG += h_cv[8];
283 }
284 }
285 GRID_LOOP_STOP(lfc, -1, 0);
286 c_M += nM;
287 a_G += nG;
288 }
289 }
290 else {
291 const double complex* c_M = (const double complex*)PyArray_DATA(c_xM_obj);
292 double complex* a_G = (double complex*)PyArray_DATA(a_xG_obj);
293 for (int x = 0; x < nx; x++) {
294 GRID_LOOP_START(lfc, q, 0) {
295
296 // In one grid loop iteration, only i2 changes.
297 int i2 = Ga % n_c[2] + beg_c[2];
298 int i1 = (Ga / n_c[2]) % n_c[1] + beg_c[1];
299 int i0 = Ga / (n_c[2] * n_c[1]) + beg_c[0];
300 // Grid point position
301 double xG = h_cv[0] * i0 + h_cv[3] * i1 + h_cv[6] * i2;
302 double yG = h_cv[1] * i0 + h_cv[4] * i1 + h_cv[7] * i2;
303 double zG = h_cv[2] * i0 + h_cv[5] * i1 + h_cv[8] * i2;
304 // Loop over grid points in current stride
305 for (int G = Ga; G < Gb; G++) {
306 // Loop over volumes at current grid point
307 for (int i = 0; i < ni; i++) {
308 // Phase of volume
309 double complex conjphase = conj(phase_i[i]);
310 LFVolume* vol = volume_i[i];
311 int M = vol->M;
312 // Check that the volume belongs to the atom in consideration later
313 int W = vol->W;
314 int nm = vol->nm;
315 int l = (nm - 1) / 2;
316
317 const bmgsspline* spline = (const bmgsspline*) \
318 &((const SplineObject*)PyList_GetItem(spline_M_obj, M))->spline;
319
320 double x = xG - pos_Wc[W][0];
321 double y = yG - pos_Wc[W][1];
322 double z = zG - pos_Wc[W][2];
323 double R_c[] = {x, y, z};
324 double r2 = x * x + y * y + z * z;
325 double r = sqrt(r2);
326 double f;
327 double dfdr;
328
329 bmgs_get_value_and_derivative(spline, r, &f, &dfdr);
330
331 // First contribution: f * d(r^l * Y)/dv
332 double fdrlYdx_m[nm];
333 if (v == 0)
334 spherical_harmonics_derivative_x(l, f, x, y, z, r2, fdrlYdx_m);
335 else if (v == 1)
336 spherical_harmonics_derivative_y(l, f, x, y, z, r2, fdrlYdx_m);
337 else
338 spherical_harmonics_derivative_z(l, f, x, y, z, r2, fdrlYdx_m);
339
340 for (int m = 0; m < nm; m++)
341 a_G[G] += fdrlYdx_m[m] * c_M[M + m] * conjphase;
342
343 // Second contribution: r^(l-1) * Y * df/dr * R_v
344 if (r > 1e-15) {
345 double rlm1Ydfdr_m[nm]; // r^(l-1) * Y * df/dr
346 double rm1dfdr = 1. / r * dfdr;
347 spherical_harmonics(l, rm1dfdr, x, y, z, r2, rlm1Ydfdr_m);
348 for (int m = 0; m < nm; m++)
349 a_G[G] += rlm1Ydfdr_m[m] * R_c[v] * c_M[M + m] * conjphase;
350 }
351 }
352 // Update coordinates of current grid point
353 xG += h_cv[6];
354 yG += h_cv[7];
355 zG += h_cv[8];
356 }
357 }
358 GRID_LOOP_STOP(lfc, q, 0);
359 c_M += nM;
360 a_G += nG;
361 }
362 }
363 Py_RETURN_NONE;
364 }
365