1 #include <Python.h>
2 #define PY_ARRAY_UNIQUE_SYMBOL GPAW_ARRAY_API
3 #define NO_IMPORT_ARRAY
4 #include <numpy/arrayobject.h>
5 #include "extensions.h"
6
7 #define VELOCITY_VERLET_ADJUST_POSITION_TOL 1e-13
8 #define VELOCITY_VERLET_ADJUST_VELOCITY_TOL 1e-13
9 #define COULOMB_CONSTANT 14.399645351950548
10 // #define FF_DEBUG_M 1
11
12 #ifdef __clang__
13 #define INLINE static
14 #else
15 #define INLINE static inline
16 #endif
17
18 // Index notation
19 // v is coordinate (x,y or z)
20 // x is distance index (1-2, 2-3, or 3-1)
21 // n is water index
22 // i is atom index (O, H1 or H2)
23
vec3_sub(double * target_v,double * a_v,double * b_v)24 INLINE void vec3_sub(double* target_v, double* a_v, double* b_v)
25 {
26 target_v[0] = a_v[0] - b_v[0];
27 target_v[1] = a_v[1] - b_v[1];
28 target_v[2] = a_v[2] - b_v[2];
29 }
30
vec3_submin(double * target_v,double * a_v,double * b_v,unsigned char * pbc,double * celldiag)31 INLINE void vec3_submin(double* target_v, double* a_v, double* b_v, unsigned char* pbc, double* celldiag)
32 {
33 for (unsigned int v=0; v<3; v++)
34 if (pbc[v]) {
35 double s = a_v[v] - b_v[v];
36 s -= round(s / celldiag[v])*celldiag[v];
37 target_v[v] = s;
38 } else {
39 target_v[v] = a_v[v] - b_v[v];
40 }
41 }
42
vec9_diffs(double * target_xv,double * R_iv)43 INLINE void vec9_diffs(double* target_xv, double* R_iv)
44 {
45 // R1-R2
46 target_xv[0] = R_iv[0*3 + 0] - R_iv[1*3 + 0]; // dOH1_x
47 target_xv[1] = R_iv[0*3 + 1] - R_iv[1*3 + 1]; // dOH1_y
48 target_xv[2] = R_iv[0*3 + 2] - R_iv[1*3 + 2]; // dOH1_z
49 // R2-R3
50 target_xv[3] = R_iv[1*3 + 0] - R_iv[2*3 + 0]; // dOH2_x
51 target_xv[4] = R_iv[1*3 + 1] - R_iv[2*3 + 1]; // dOH2_y
52 target_xv[5] = R_iv[1*3 + 2] - R_iv[2*3 + 2]; // dOH2_z
53 // R3-R1
54 target_xv[6] = R_iv[2*3 + 0] - R_iv[0*3 + 0]; // dH1H2_x
55 target_xv[7] = R_iv[2*3 + 1] - R_iv[0*3 + 1]; // dH1H2_y
56 target_xv[8] = R_iv[2*3 + 2] - R_iv[0*3 + 2]; // dH1H2_z
57 }
58
vec9_massdiffs(double * target_xv,double * im_i,double * P_iv)59 INLINE void vec9_massdiffs(double* target_xv, double* im_i, double* P_iv)
60 {
61 // R1-R2
62 target_xv[0] = im_i[0] * P_iv[0*3 + 0] - im_i[1] * P_iv[1*3 + 0]; // dOH1_x
63 target_xv[1] = im_i[0] * P_iv[0*3 + 1] - im_i[1] * P_iv[1*3 + 1]; // dOH1_y
64 target_xv[2] = im_i[0] * P_iv[0*3 + 2] - im_i[1] * P_iv[1*3 + 2]; // dOH1_z
65 // R2-R3
66 target_xv[3] = im_i[1] * P_iv[1*3 + 0] - im_i[2] * P_iv[2*3 + 0]; // dOH2_x
67 target_xv[4] = im_i[1] * P_iv[1*3 + 1] - im_i[2] * P_iv[2*3 + 1]; // dOH2_y
68 target_xv[5] = im_i[1] * P_iv[1*3 + 2] - im_i[2] * P_iv[2*3 + 2]; // dOH2_z
69 // R3-R1
70 target_xv[6] = im_i[2] * P_iv[2*3 + 0] - im_i[0] * P_iv[0*3 + 0]; // dH1H2_x
71 target_xv[7] = im_i[2] * P_iv[2*3 + 1] - im_i[0] * P_iv[0*3 + 1]; // dH1H2_y
72 target_xv[8] = im_i[2] * P_iv[2*3 + 2] - im_i[0] * P_iv[0*3 + 2]; // dH1H2_z
73 }
74
vec9_dot(double * lambda_x,double * newd_xv,double * d_xv)75 INLINE void vec9_dot(double* lambda_x, double* newd_xv, double* d_xv)
76 {
77 lambda_x[0] = newd_xv[0*3 + 0] * d_xv[0*3 + 0] +
78 newd_xv[0*3 + 1] * d_xv[0*3 + 1] +
79 newd_xv[0*3 + 2] * d_xv[0*3 + 2];
80 lambda_x[1] = newd_xv[1*3 + 0] * d_xv[1*3 + 0] +
81 newd_xv[1*3 + 1] * d_xv[1*3 + 1] +
82 newd_xv[1*3 + 2] * d_xv[1*3 + 2];
83 lambda_x[2] = newd_xv[2*3 + 0] * d_xv[2*3 + 0] +
84 newd_xv[2*3 + 1] * d_xv[2*3 + 1] +
85 newd_xv[2*3 + 2] * d_xv[2*3 + 2];
86 }
87
vec3_imul(double * target_x,double * a_x)88 INLINE void vec3_imul(double* target_x, double* a_x)
89 {
90 target_x[0]*= a_x[0];
91 target_x[1]*= a_x[1];
92 target_x[2]*= a_x[2];
93 }
94
vec3_div(double * target_x,double * a_x,double * b_x)95 INLINE void vec3_div(double* target_x, double* a_x, double* b_x)
96 {
97 target_x[0]= a_x[0] / b_x[0];
98 target_x[1]= a_x[1] / b_x[1];
99 target_x[2]= a_x[2] / b_x[2];
100 }
101
vec3_axpy(double * target_v,double coeff,double * a_v)102 INLINE void vec3_axpy(double* target_v, double coeff, double* a_v)
103 {
104 target_v[0] += coeff * a_v[0];
105 target_v[1] += coeff * a_v[1];
106 target_v[2] += coeff * a_v[2];
107 }
108
sqr(double a)109 INLINE double sqr(double a)
110 {
111 return a*a;
112 }
113
vec9_sqrsum(double * target_x,double * R_xv)114 INLINE void vec9_sqrsum(double* target_x, double* R_xv)
115 {
116 target_x[0] = sqr(R_xv[0*3 + 0]) + sqr(R_xv[0*3 + 1]) + sqr(R_xv[0*3 + 2]);
117 target_x[1] = sqr(R_xv[1*3 + 0]) + sqr(R_xv[1*3 + 1]) + sqr(R_xv[1*3 + 2]);
118 target_x[2] = sqr(R_xv[2*3 + 0]) + sqr(R_xv[2*3 + 1]) + sqr(R_xv[2*3 + 2]);
119 }
120
vec3_sqrsum(double * R_v)121 INLINE double vec3_sqrsum(double* R_v)
122 {
123 return sqr(R_v[0]) + sqr(R_v[1]) + sqr(R_v[2]);
124 }
125
vec9_print(char * title,double * target_iv)126 void vec9_print(char* title, double* target_iv)
127 {
128 printf("%s:\n", title);
129 printf("%20.15f %20.15f %20.15f\n", target_iv[0], target_iv[1], target_iv[2]);
130 printf("%20.15f %20.15f %20.15f\n", target_iv[3], target_iv[4], target_iv[5]);
131 printf("%20.15f %20.15f %20.15f\n", target_iv[6], target_iv[7], target_iv[8]);
132 }
133
vec3_print(char * title,double * target_v)134 void vec3_print(char* title, double* target_v)
135 {
136 printf("%s:\n", title);
137 printf("%20.15f %20.15f %20.15f\n", target_v[0], target_v[1], target_v[2]);
138 }
139
coulomb(double Za,double Zb,double * Ra_v,double * Rb_v,double * Fa_v,double * Fb_v,unsigned char * pbc,double * celldiag)140 INLINE double coulomb(double Za, double Zb, double* Ra_v, double* Rb_v, double* Fa_v, double* Fb_v, unsigned char* pbc, double* celldiag)
141 {
142 double d_v[3];
143 vec3_submin(d_v, Ra_v, Rb_v, pbc, celldiag);
144 double r2 = vec3_sqrsum(d_v);
145 double r = sqrt(r2);
146 double E = COULOMB_CONSTANT * Za * Zb / r;
147 double str = E / r2;
148 vec3_axpy(Fa_v, str, d_v);
149 vec3_axpy(Fb_v, -str, d_v);
150 return E;
151 }
152
coulomb_cutoff(double Za,double Zb,double t,double * Ra_v,double * Rb_v,double * Fa_v,double * Fb_v,unsigned char * pbc,double * celldiag)153 INLINE double coulomb_cutoff(double Za, double Zb, double t, double* Ra_v, double* Rb_v, double* Fa_v, double* Fb_v, unsigned char* pbc, double* celldiag)
154 {
155 double d_v[3];
156 vec3_submin(d_v, Ra_v, Rb_v, pbc, celldiag);
157 double r2 = vec3_sqrsum(d_v);
158 double r = sqrt(r2);
159 double E = COULOMB_CONSTANT * Za * Zb / r;
160 double str = E / r2 * t;
161 vec3_axpy(Fa_v, str, d_v);
162 vec3_axpy(Fb_v, -str, d_v);
163 return E; // The energy is not correct energy, it is multiplied later with t
164 }
165
LJ(double A,double B,double * Ra_v,double * Rb_v,double * Fa_v,double * Fb_v,unsigned char * pbc,double * celldiag)166 INLINE double LJ(double A, double B, double* Ra_v, double* Rb_v, double* Fa_v, double* Fb_v, unsigned char* pbc, double* celldiag)
167 {
168 double d_v[3];
169 vec3_submin(d_v, Ra_v, Rb_v, pbc, celldiag);
170 double r2 = vec3_sqrsum(d_v);
171 double r4 = r2*r2;
172 double r6 = r4*r2;
173 double r8 = r4*r4;
174 double r12 = r8*r4;
175 double r14 = r12*r2;
176 double str = 12*A/r14+6*B/r8;
177 vec3_axpy(Fa_v, str, d_v);
178 vec3_axpy(Fb_v, -str, d_v);
179 return A / r12 + B/r6;
180 }
181
pair_interaction_cutoff(double A,double B,double cutoff,double width,double * Z_i,double * Ra_iv,double * Rb_iv,double * Fa_iv,double * Fb_iv,unsigned char * pbc,double * celldiag)182 INLINE double pair_interaction_cutoff(double A, double B, double cutoff, double width, double* Z_i, double* Ra_iv, double* Rb_iv, double* Fa_iv, double* Fb_iv, unsigned char* pbc, double* celldiag)
183 {
184 double d_v[3];
185 vec3_submin(d_v, Ra_iv + 0*3, Rb_iv + 0*3, pbc, celldiag);
186 double r2 = vec3_sqrsum(d_v);
187 double t;
188 double dtdr;
189 double r;
190 if (r2 > (cutoff*cutoff)) { return 0.0; }
191 if (r2 < (cutoff-width)*(cutoff-width)) {
192 t = 1;
193 dtdr = 0;
194 r = 1; // XXX Dummy value, this is not used because dtdr = 0
195 } else {
196 r = sqrt(r2);
197 double y = (r - cutoff + width) / width;
198 t = 1.0-y*y * (3.0 - 2.0 * y);
199 dtdr = -6.0 / width * y * (1.0 - y);
200 }
201
202 double r4 = r2*r2;
203 double r6 = r4*r2;
204 double r8 = r4*r4;
205 double r12 = r8*r4;
206 double r14 = r12*r2;
207 double E = (A / r12 + B/r6);
208 double str = (12*A/r14+6*B/r8)* t - E*dtdr / r;
209 E *= t;
210 vec3_axpy(Fa_iv + 0*3, str, d_v);
211 vec3_axpy(Fb_iv + 0*3, -str, d_v);
212
213 for (unsigned int i1 = 0; i1 < 3; i1++)
214 for (unsigned int i2 = 0; i2 < 3; i2++) {
215 double Ept = coulomb_cutoff(Z_i[i1], Z_i[i2], t, Ra_iv+i1*3, Rb_iv+i2*3, Fa_iv+i1*3, Fb_iv+i2*3, pbc, celldiag);
216 E+= Ept * t;
217 double str = Ept / r * dtdr;
218 vec3_axpy(Fa_iv + 0*3, -str, d_v); // Note that these are added to O, even for H-H interaction.
219 vec3_axpy(Fb_iv + 0*3, str, d_v);
220 }
221
222 return E;
223
224 }
225
adjust_positions(PyObject * self,PyObject * args)226 PyObject* adjust_positions(PyObject *self, PyObject *args)
227 {
228 PyArrayObject* arraylen_x = 0; // Input: the 3 constraint lengths
229 PyArrayObject* arraymass_i = 0; // Input: the 3 masses
230 PyArrayObject* arrayR_niv = 0; // Output: Adjusted positions will be written here.
231 PyArrayObject* arraynewR_niv = 0; // Input: gives the positions to be adjusted.
232
233 if (!PyArg_ParseTuple(args, "OOOO", &arraylen_x, &arraymass_i,
234 &arrayR_niv, &arraynewR_niv))
235 {
236 return NULL;
237 }
238
239 unsigned int NA = PyArray_DIM(arrayR_niv, 0);
240
241 if (NA % 3 != 0)
242 {
243 PyErr_SetString(PyExc_TypeError, "Number of atoms not divisible with 3.");
244 return NULL;
245 }
246
247 unsigned int NW = NA / 3;
248
249 if (!(PyArray_NDIM(arraymass_i) == 1 &&
250 PyArray_DIM(arraymass_i,0) == 3))
251 {
252 PyErr_SetString(PyExc_TypeError, "mass_i should be array with length 3.");
253 return NULL;
254 }
255
256 if (!(PyArray_NDIM(arraylen_x) == 1 &&
257 PyArray_DIM(arraylen_x,0) == 3))
258 {
259 PyErr_SetString(PyExc_TypeError, "len_x should be array with length 3.");
260 return NULL;
261 }
262
263 double* len_x = DOUBLEP(arraylen_x);
264 double* mass_i = DOUBLEP(arraymass_i);
265
266 double len2_x[3];
267 len2_x[0] = sqr(len_x[0]);
268 len2_x[1] = sqr(len_x[1]);
269 len2_x[2] = sqr(len_x[2]);
270
271 double mu_x[3]; // Reduced masses of pairs
272 mu_x[0] = 1.0 / ( (1.0/mass_i[0]) + (1.0/mass_i[1]) );
273 mu_x[1] = 1.0 / ( (1.0/mass_i[1]) + (1.0/mass_i[2]) );
274 mu_x[2] = 1.0 / ( (1.0/mass_i[2]) + (1.0/mass_i[0]) );
275
276 double invm_i[3]; // Inverse masses of atoms
277 invm_i[0] = 0.5 / mass_i[0]; // Note: /2 is embedded here
278 invm_i[1] = 0.5 / mass_i[1];
279 invm_i[2] = 0.5 / mass_i[2];
280
281 double* R_niv = DOUBLEP(arrayR_niv);
282 double* newR_niv = DOUBLEP(arraynewR_niv);
283
284 for (unsigned int n=0; n<NW; n++) // For each water molecule
285 {
286 #ifdef FF_DEBUG
287 vec9_print("R_iv", newR_niv+n*9);
288 #endif
289 double d_xv[9];
290 // 1. Calculate (R1-R2), (R2-R3) and (R3-R1)
291 vec9_diffs(d_xv, R_niv+n*9);
292 #ifdef FF_DEBUG
293 vec9_print("d_xv", d_xv);
294 #endif
295
296 // 2. Calculate (R1-R2)^2, (R2-R3)^2, (R3-R1)^2
297 // double R_x[3];
298 // vec9_sqrsum(R_x, d_xv);
299
300 int iteration = 0;
301 while (1) {
302 #ifdef FF_DEBUG
303 printf("Iteration: %d\n", iteration);
304 vec9_print("newR_iv", newR_niv+n*9);
305 #endif
306 double newd_xv[9];
307 // 1. Calculate (R1'-R2'), (R2'-R3') and (R3'-R1')
308 vec9_diffs(newd_xv, newR_niv+n*9);
309 #ifdef FF_DEBUG
310 vec9_print("newd_xv", newd_xv);
311 #endif
312
313 // 2. Calculate (R1'-R2')^2, (R2'-R3')^2, (R3'-R1')^2
314 double newR_x[3];
315 vec9_sqrsum(newR_x, newd_xv);
316 #ifdef FF_DEBUG
317 vec3_print("newR_x", newR_x);
318 #endif
319
320 // 3. Calculate f12 = (R1'-R2')^2-D12^2, f23 and f31
321 double f_x[3];
322 vec3_sub(f_x, newR_x, len2_x);
323 #ifdef FF_DEBUG
324 vec3_print("f_x", f_x);
325 #endif
326
327 if (iteration++ > 1000) {
328 printf("Warning: Adjust positions did not converge.\n");
329 break;
330 }
331 if ((fabs(f_x[0]) < VELOCITY_VERLET_ADJUST_POSITION_TOL) &&
332 (fabs(f_x[1]) < VELOCITY_VERLET_ADJUST_POSITION_TOL) &&
333 (fabs(f_x[2]) < VELOCITY_VERLET_ADJUST_POSITION_TOL)) {
334 break;
335 }
336
337 double lambda_x[3];
338 double denom_x[3];
339 // Calculate lambdas
340 vec9_dot(denom_x, newd_xv, d_xv); // (R1-R2) . (R1'-R2') etc.
341 vec3_div(lambda_x, f_x, denom_x); // f_12 / ( (R1-R2) . (R1'-R2') etc.
342 vec3_imul(lambda_x, mu_x); // lambda /= (m_1^-1 + m_2^-1) etc.
343 #ifdef FF_DEBUG
344 vec3_print("lambda_x", lambda_x);
345 #endif
346
347 // Update newR's
348 // newR_1 += -lambda12 * (R1-R2)
349 vec3_axpy(newR_niv + n*9 + 0*3, -lambda_x[0] * invm_i[0], d_xv + 0*3);
350 // newR_1 += lambda31 * (R3-R1)
351 vec3_axpy(newR_niv + n*9 + 0*3, lambda_x[2] * invm_i[0], d_xv + 2*3);
352
353 // newR_2 += +lambda12 * (R1-R2)
354 vec3_axpy(newR_niv + n*9 + 1*3, +lambda_x[0] * invm_i[1], d_xv + 0*3);
355 // newR_2 += -lambda23 * (R2-R3)
356 vec3_axpy(newR_niv + n*9 + 1*3, -lambda_x[1] * invm_i[1], d_xv + 1*3);
357
358 // newR_3 += +lambda23 * (R2-R3)
359 vec3_axpy(newR_niv + n*9 + 2*3, +lambda_x[1] * invm_i[2], d_xv + 1*3);
360 // newR_3 += -lambda31 * (R3-R1)
361 vec3_axpy(newR_niv + n*9 + 2*3, -lambda_x[2] * invm_i[2], d_xv + 2*3);
362 }
363 }
364 Py_RETURN_NONE;
365 }
366
adjust_momenta(PyObject * self,PyObject * args)367 PyObject* adjust_momenta(PyObject *self, PyObject *args)
368 {
369 PyArrayObject* arraymass_i = 0; // Input: the 3 masses
370 PyArrayObject* arrayR_niv = 0; //
371 PyArrayObject* arraynewP_niv = 0; //
372
373 if (!PyArg_ParseTuple(args, "OOO", &arraymass_i, &arrayR_niv,
374 &arraynewP_niv))
375 {
376 return NULL;
377 }
378
379 unsigned int NA = PyArray_DIM(arrayR_niv, 0);
380
381 if (NA % 3 != 0)
382 {
383 PyErr_SetString(PyExc_TypeError, "Number of atoms not divisible with 3.");
384 return NULL;
385 }
386
387 unsigned int NW = NA / 3;
388
389 if (!(PyArray_NDIM(arraymass_i) == 1 &&
390 PyArray_DIM(arraymass_i,0) == 3))
391 {
392 PyErr_SetString(PyExc_TypeError, "mass_i should be array with length 3.");
393 return NULL;
394 }
395
396 double* mass_i = DOUBLEP(arraymass_i);
397
398 double invm_i[3]; // Inverse masses of atoms
399 invm_i[0] = 1.0 / mass_i[0];
400 invm_i[1] = 1.0 / mass_i[1];
401 invm_i[2] = 1.0 / mass_i[2];
402
403 double mu_x[3]; // Reduced masses of pairs
404 mu_x[0] = 1.0 / ( (1.0/mass_i[0]) + (1.0/mass_i[1]) );
405 mu_x[1] = 1.0 / ( (1.0/mass_i[1]) + (1.0/mass_i[2]) );
406 mu_x[2] = 1.0 / ( (1.0/mass_i[2]) + (1.0/mass_i[0]) );
407
408 double* R_niv = DOUBLEP(arrayR_niv);
409 double* newP_niv = DOUBLEP(arraynewP_niv);
410
411 for (unsigned int n=0; n<NW; n++) // For each water molecule
412 {
413 #ifdef FF_DEBUG_M
414 vec9_print("R_iv", newP_niv+n*9);
415 #endif
416 double d_xv[9];
417 // 1. Calculate (R1-R2), (R2-R3) and (R3-R1)
418 vec9_diffs(d_xv, R_niv+n*9);
419 #ifdef FF_DEBUG_M
420 vec9_print("d_xv", d_xv);
421 #endif
422
423 int iteration = 0;
424 while (1) {
425 #ifdef FF_DEBUG_M
426 printf("Iteration: %d\n", iteration);
427 vec9_print("newP_iv", newP_niv+n*9);
428 #endif
429 double newd_xv[9];
430 // 1. Calculate (m_1^-1 P1'-m_2^-1 P2'), etc.
431 vec9_massdiffs(newd_xv, invm_i, newP_niv+n*9);
432
433 #ifdef FF_DEBUG_M
434 vec9_print("(m_1^-1 P1'-m_2^-1 P2')", newd_xv);
435 #endif
436
437 // 3. Calculate g12 = (R1-R2) . (m_1^-1 P1'-m_2^-1 P2') etc.
438 double g_x[3];
439 vec9_dot(g_x, newd_xv, d_xv);
440
441 #ifdef FF_DEBUG_M
442 vec3_print("g_x", g_x);
443 #endif
444
445 if (iteration++ > 1000) {
446 printf("Warning: Adjust velocities did not converge.\n");
447 break;
448 }
449
450 if ((fabs(g_x[0]) < VELOCITY_VERLET_ADJUST_VELOCITY_TOL) &&
451 (fabs(g_x[1]) < VELOCITY_VERLET_ADJUST_VELOCITY_TOL) &&
452 (fabs(g_x[2]) < VELOCITY_VERLET_ADJUST_VELOCITY_TOL)) {
453 break;
454 }
455
456 double lambda_x[3];
457 double denom_x[3];
458
459 vec9_dot(denom_x, d_xv, d_xv); // (R1-R2)^2 etc.
460 vec3_div(lambda_x, g_x, denom_x); // g_12 / (R1-R2)^2
461 vec3_imul(lambda_x, mu_x); // lambda /= (m_1^-1 + m_2^-1) etc.
462
463 #ifdef FF_DEBUG_M
464 vec3_print("lambda_x", lambda_x);
465 #endif
466
467 // Update newR's
468 // newR_1 += -lambda12 * (R1-R2)
469 vec3_axpy(newP_niv + n*9 + 0*3, -lambda_x[0], d_xv + 0*3);
470 // newR_1 += lambda31 * (R3-R1)
471 vec3_axpy(newP_niv + n*9 + 0*3, lambda_x[2], d_xv + 2*3);
472
473 // newR_2 += +lambda12 * (R1-R2)
474 vec3_axpy(newP_niv + n*9 + 1*3, +lambda_x[0], d_xv + 0*3);
475 // newR_2 += -lambda23 * (R2-R3)
476 vec3_axpy(newP_niv + n*9 + 1*3, -lambda_x[1], d_xv + 1*3);
477
478 // newR_3 += +lambda23 * (R2-R3)
479 vec3_axpy(newP_niv + n*9 + 2*3, +lambda_x[1], d_xv + 1*3);
480 // newR_3 += -lambda31 * (R3-R1)
481 vec3_axpy(newP_niv + n*9 + 2*3, -lambda_x[2], d_xv + 2*3);
482 }
483 }
484 Py_RETURN_NONE;
485 }
486
calculate_forces_H2O(PyObject * self,PyObject * args)487 PyObject* calculate_forces_H2O(PyObject *self, PyObject *args)
488 {
489 double A;
490 double B;
491 PyArrayObject* arraypbc = 0;
492 PyArrayObject* arrayZ_i = 0;
493 PyArrayObject* arrayR_niv = 0;
494 PyArrayObject* arrayF_niv = 0;
495 PyArrayObject* arraycell = 0;
496
497 double cutoff=0;
498 double width=0;
499
500 if (!PyArg_ParseTuple(args, "OOddddOOO", &arraypbc, &arraycell, &A, &B, &cutoff, &width, &arrayZ_i, &arrayR_niv, &arrayF_niv))
501 {
502 return NULL;
503 }
504
505 if (!(PyArray_NDIM(arraypbc) == 1 &&
506 PyArray_DIM(arraypbc,0) == 3)) {
507 PyErr_SetString(PyExc_TypeError, "pbc should be array with length 3.");
508 return NULL;
509 }
510
511 if (!(PyArray_NDIM(arrayZ_i) == 1 &&
512 PyArray_DIM(arrayZ_i,0) == 3)) {
513 PyErr_SetString(PyExc_TypeError, "Z_i should be array with length 3.");
514 return NULL;
515 }
516
517 if (!(PyArray_NDIM(arraycell) == 2 &&
518 (PyArray_DIM(arraycell,0) == 3) &&
519 (PyArray_DIM(arraycell,1) == 3))) {
520 PyErr_SetString(PyExc_TypeError, "Cell should be array with size 3x3.");
521 return NULL;
522 }
523
524 double* cell_vc = DOUBLEP(arraycell);
525 //printf("Cell\n");
526 for (unsigned int v1=0; v1<3; v1++)
527 for (unsigned int v2=0; v2<3; v2++) {
528 //printf("%20.15f \n", cell_vc[v1+v2*3]);
529 if (v1 != v2)
530 if (fabs(cell_vc[v1+v2*3]) > 1e-10) {
531 PyErr_SetString(PyExc_TypeError, "Cell array should be diagonal.");
532 return NULL;
533 }
534 }
535
536 double celldiag[3];
537 celldiag[0] = cell_vc[0];
538 celldiag[1] = cell_vc[1+1*3];
539 celldiag[2] = cell_vc[2+2*3];
540
541 unsigned int NA = PyArray_DIM(arrayR_niv, 0);
542
543 if (NA % 3 != 0)
544 {
545 PyErr_SetString(PyExc_TypeError, "Number of atoms not divisible with 3.");
546 return NULL;
547 }
548
549 unsigned char *pbc = (unsigned char*) PyArray_DATA(arraypbc);
550 //printf("Boundary conditions %u %u %u \n", pbc[0], pbc[1], pbc[2]);
551
552 unsigned int NW = NA / 3;
553
554 double* Z_i = DOUBLEP(arrayZ_i);
555 double* R_niv = DOUBLEP(arrayR_niv);
556 double* F_niv = DOUBLEP(arrayF_niv);
557
558 double E = 0.0;
559 if (cutoff <= 0.0) { // No cutoff
560 // Double loop of atoms
561 for (unsigned int n1 = 0; n1 < NW; n1++)
562 for (unsigned int n2 = n1+1; n2 < NW; n2++) {
563 E += LJ(A, B, R_niv+n1*9, R_niv+n2*9, F_niv+n1*9, F_niv+n2*9, pbc, celldiag);
564 for (unsigned int i1 = 0; i1 < 3; i1++)
565 for (unsigned int i2 = 0; i2 < 3; i2++)
566 E+= coulomb(Z_i[i1], Z_i[i2], R_niv+n1*9+i1*3, R_niv+n2*9+i2*3, F_niv+n1*9+i1*3, F_niv+n2*9+i2*3, pbc, celldiag);
567 }
568 } else { // With cutoff
569 // Double loop of atoms
570 for (unsigned int n1 = 0; n1 < NW; n1++)
571 for (unsigned int n2 = n1+1; n2 < NW; n2++) {
572 E += pair_interaction_cutoff(A, B, cutoff, width, Z_i, R_niv+n1*9, R_niv+n2*9, F_niv+n1*9, F_niv+n2*9, pbc, celldiag);
573 }
574 }
575
576 return Py_BuildValue("d", E);
577 }
578
579