1 """This modules implements the mathematics-related classes, functions and procedures.
2
3 Classes:
4 * QVectors: the class that actually performs the q vectors generation.
5
6 Functions:
7 * differentiate : performs a numerical differentiation of 1D Numeric array
8 * correlation : performs the numerical correlation between two 1D Numeric arrays
9 * convolution : performs the numerical convolution between two 1D Numeric arrays
10 * FFT : performs the FFT of a 1D Numeric array
11 * invFFT : performs the inverse FFT of a 1D Numeric array
12 * gaussianWindow : performs a Gaussian smoothing of 1D Numeric array.
13 * factorial : computes factorial (n) where n is an integer.
14 * basisVectors : computes the basis vectors of the simulation cell from a set of values defining its geometry (3 distances and 3 angles).
15 * randomPointInCircle: returns a vector within a circle of radius |r| and orthogonal to a given direction.
16 * randomDirection2D : returns a normalized vector generated from a unit circle orthogonal to a given direction.
17 * randomPlane2D : generates a normalized random q-vector on a plane defined by vect1, vect2.
18 * qVectorGenerator : sets up and returns a set of q vectors generated from different user-defined parameters.
19 * sphericalHarmonics : calculates the spherical functions Y from a set of j, m, n Wigner indexes
20 * preparePP : sets up the calculation of spherical harmonics.
21 """
22
23
24 from random import uniform
25 import sys
26
27
28 from Scientific import N
29 from Scientific.FFT import fft, inverse_fft
30 from Scientific.Geometry import Tensor, Vector
31
32
33 from MMTK.Random import randomDirection
34
35
36 from nMOLDYN.Core.Error import Error
37 from nMOLDYN.Core.Logger import LogMessage
38
39
40
41
42
43 """a2 = array used to perform order 2 numerical differentiation scheme."""
44 a2 = N.array([[ -3., 4., -1.],
45 [ -1., 0., 1.],
46 [ 1., -4., 3.]])
47
48 """a3 = array used to perform order 3 numerical differentiation scheme."""
49 a3 = N.array([[ -11., 18., -9., 2.],
50 [ -2., -3., 6., -1.],
51 [ 1., -6., 3., 2.],
52 [ -2., 9., -18., 11.]])
53
54 """a4 = array used to perform order 4 numerical differentiation scheme."""
55 a4 = N.array([[ -50., 96., -72., 32., -6.],
56 [ -6., -20., 36., -12., 2.],
57 [ 2., -16., 0., 16., -2.],
58 [ -2., 12., -36., 20., 6.],
59 [ 6., -32., 72., -96., 50.]])
60
61 """a5 = N.array used to perform order 5 numerical differentiation scheme."""
62 a5 = N.array([[ -274., 600., -600., 400., -150., 24.],
63 [ -24., -130., 240., -120., 40., -6.],
64 [ 6., -60., -40., 120., -30., 4.],
65 [ -4., 30., -120., 40., 60., -6.],
66 [ 6., -40., 120., -240., 130., 24.],
67 [ -24., 150., -400., 600., -600., 274.]])
68
70 """Returns the numerical derivative of order |order| of the signal |inputSeries| using the
71 differentiation step |dx|.
72
73 @param inputSeries: the signal to differentiate.
74 @index inputSeries: NumPy array.
75
76 @param order: an integer in [1,5] specifying the numerical differentiation order.
77 @index order: integer
78
79 @param dx: a float specifying the differentiation step. Assumed to be constant over all the spectrum.
80 @index dx: float
81
82 @return: the differentiated signal.
83 @rtype: NumPy array
84
85 @see: M. Abramowitz, I.A. Stegun; 'Handbook of mathematical functions', Dover, New-York, 1972 p.914.
86 """
87
88 if len(inputSeries) == 0:
89 raise Error('The inputSeries is empty.')
90
91 if dx <= 0.0:
92 raise Error('The time step for derivation must strictly positive.')
93
94
95 outputSeries = N.zeros(len(inputSeries),'d')
96
97 if order == 1:
98
99
100
101
102
103
104 fact = 1.0/(2.0*dx)
105
106
107 gj = inputSeries[1:] - inputSeries[:-1]
108
109
110 outputSeries[0] = N.add.reduce(a2[0,:]*inputSeries[:3])*fact
111 outputSeries[-1] = N.add.reduce(a2[2,:]*inputSeries[-3:])*fact
112
113
114
115
116 outputSeries[1:-1] = (gj[1:]+gj[:-1])*fact
117
118 elif order == 2:
119
120
121
122 fact = 1.0/(2.0*dx)
123
124
125 outputSeries[0] = N.add.reduce(a2[0,:]*inputSeries[:3])*fact
126 outputSeries[-1] = N.add.reduce(a2[2,:]*inputSeries[-3:])*fact
127
128
129 gj = N.zeros((len(inputSeries)-2,3),'d')
130 gj[:,0] = a2[1,0]*inputSeries[:-2]
131 gj[:,1] = a2[1,1]*inputSeries[1:-1]
132 gj[:,2] = a2[1,2]*inputSeries[2:]
133 outputSeries[1:-1] = N.add.reduce(gj,-1)*fact
134
135 elif order == 3:
136
137
138
139 fact = 1./(6.*dx)
140
141
142 outputSeries[0] = N.add.reduce(a3[0,:]*inputSeries[:4])*fact
143 outputSeries[1] = N.add.reduce(a3[1,:]*inputSeries[:4])*fact
144 outputSeries[-1] = N.add.reduce(a3[3,:]*inputSeries[-4:])*fact
145
146
147 gj = N.zeros((len(inputSeries)-3,4),'d')
148 gj[:,0] = a3[2,0]*inputSeries[:-3]
149 gj[:,1] = a3[2,1]*inputSeries[1:-2]
150 gj[:,2] = a3[2,2]*inputSeries[2:-1]
151 gj[:,3] = a3[2,3]*inputSeries[3:]
152 outputSeries[2:-1] = N.add.reduce(gj,-1)*fact
153
154 elif order == 4:
155
156
157
158 fact = 1./(24.*dx)
159
160
161 outputSeries[0] = N.add.reduce(a4[0,:]*inputSeries[:5])*fact
162 outputSeries[1] = N.add.reduce(a4[1,:]*inputSeries[:5])*fact
163 outputSeries[-2] = N.add.reduce(a4[3,:]*inputSeries[-5:])*fact
164 outputSeries[-1] = N.add.reduce(a4[4,:]*inputSeries[-5:])*fact
165
166
167 gj = N.zeros((len(inputSeries)-4,5),'d')
168 gj[:,0] = a4[2,0]*inputSeries[:-4]
169 gj[:,1] = a4[2,1]*inputSeries[1:-3]
170 gj[:,2] = a4[2,2]*inputSeries[2:-2]
171 gj[:,3] = a4[2,3]*inputSeries[3:-1]
172 gj[:,4] = a4[2,4]*inputSeries[4:]
173 outputSeries[2:-2] = N.add.reduce(gj,-1)*fact
174
175 elif order == 5:
176
177
178
179 fact = 1./(120.*dx)
180
181
182 outputSeries[0] = N.add.reduce(a5[0,:]*inputSeries[:6])*fact
183 outputSeries[1] = N.add.reduce(a5[1,:]*inputSeries[:6])*fact
184 outputSeries[2] = N.add.reduce(a5[2,:]*inputSeries[:6])*fact
185 outputSeries[-2] = N.add.reduce(a5[4,:]*inputSeries[-6:])*fact
186 outputSeries[-1] = N.add.reduce(a5[5,:]*inputSeries[-6:])*fact
187
188
189 gj = N.zeros((len(inputSeries)-5,6),'d')
190 gj[:,0] = a5[3,0]*inputSeries[:-5]
191 gj[:,1] = a5[3,1]*inputSeries[1:-4]
192 gj[:,2] = a5[3,2]*inputSeries[2:-3]
193 gj[:,3] = a5[3,3]*inputSeries[3:-2]
194 gj[:,4] = a5[3,4]*inputSeries[4:-1]
195 gj[:,5] = a5[3,5]*inputSeries[5:]
196 outputSeries[3:-2] = N.add.reduce(gj,-1)*fact
197
198 else:
199 raise Error('Unknown differentiation order.')
200
201 return outputSeries
202
203
205
206 nTimes = len(inputSeries)
207 corr = N.zeros((nTimes,), N.Float)
208 for i in range(nTimes):
209 s = 0.0
210 for j in range(nTimes-i):
211 try:
212 s += sum(inputSeries[j]*inputSeries[j+i])
213 except:
214 s += inputSeries[j]*inputSeries[j+i]
215 corr[i] = s/float(nTimes-i)
216
217 return corr
218
220 """Returns the numerical correlation between |inputSeries1| and |inputSeries2| multidimensional
221 NumPy arrays.
222
223 @param inputSeries1: the first signal
224 @type inputSeries1: NumPy array
225
226 @param inputSeries2: if not None, the second signal to correlate with |inputSeries1| otherwise the
227 correlation will be an autocorrelation.
228 @type inputSeries2: NumPy array
229
230 @return: an array (length(|inputSeries1|)) storing the result of the correlation.
231 @rtype: NumPy array
232
233 @note: if |inputSeries1| is a multidimensional array the correlation calculation is performed on
234 the first dimension.
235
236 @note: The correlation is computed using the FCA algorithm.
237 """
238
239 if len(inputSeries1) <= 0:
240 raise Error('One or both time series are empty.')
241
242
243 inputSeries1Length = len(inputSeries1)
244
245
246 extendedLength = 2*inputSeries1Length
247
248
249
250
251
252 FFTSeries1 = fft(inputSeries1,extendedLength,0)
253
254 if inputSeries2 is None:
255
256 FFTSeries2 = FFTSeries1
257 else:
258
259
260 FFTSeries2 = fft(inputSeries2,extendedLength,0)
261
262
263 FFTSeries1 = N.conjugate(FFTSeries1)*FFTSeries2
264
265
266
267 FFTSeries1 = inverse_fft(FFTSeries1,len(FFTSeries1),0)
268
269
270
271 if len(FFTSeries1.shape) == 1:
272 corr = FFTSeries1.real[:inputSeries1Length] / (inputSeries1Length-N.arange(inputSeries1Length))
273 else:
274 corr = N.add.reduce(FFTSeries1.real[:inputSeries1Length],1) / (inputSeries1Length-N.arange(inputSeries1Length))
275
276 return corr
277
279 """Returns the numerical convolution between |inputSeries1| and |inputSeries2| one-dimensional
280 NumPy arrays.
281
282 @param inputSeries1: the first signal
283 @type inputSeries1: NumPy array
284
285 @param inputSeries2: the second signal to convolve with |inputSeries1|.
286 @type inputSeries2: NumPy array
287
288 @return: an array (length(|inputSeries1|)) storing the result of the convolution.
289 @rtype: NumPy array
290
291 @note: if |inputSeries1| is a multidimensional array the convolution calculation is performed on
292 the first dimension.
293
294 @note: the convolution is computed using the convolve function of NumPy package.
295 """
296
297 if (len(inputSeries1) == 0) | (len(inputSeries2) == 0):
298 raise Error('One or both inputSeries are empty.')
299
300
301
302 return N.convolve(inputSeries1,inputSeries2)
303
304 -def FFT(inputSeries):
305 """Returns the FFT of |inputSeries| multidimensional NumPy array.
306
307 @param inputSeries: the array on which to computes the FFT.
308 @type inputSeries: NumPy array
309
310 @return: the FFT transformed array.
311 @rtype: NumPy array
312
313 @note: the FFT is computed using the fft function of Scientific.FFT package.
314 """
315
316 if len(inputSeries) == 0:
317 raise Error('The inputSeries is empty.')
318
319
320 return fft(inputSeries,len(inputSeries),0)
321
323 """Returns the inverse FFT of |inputSeries| multidimensional NumPy array.
324
325 @param inputSeries: the array on which to computes the inverse FFT.
326 @type inputSeries: NumPy array
327
328 @return: the inverse FFT transformed array.
329 @rtype: NumPy array
330
331 @note: the inverse FFT is computed using the inverse_fft function of Scientific.FFT package.
332 """
333
334 if len(inputSeries) == 0:
335 raise Error('The inputSeries is empty.')
336
337
338 return inverse_fft(inputSeries,len(inputSeries),0)
339
341 """Returns a smoothed signal using |inputSeries| input signal and a gaussian kernel of width |alpha|.
342
343 @param inputSeries: the signal to smooth.
344 @type inputSeries: NumPy array
345
346 @param alpha: a float specifying the width of the Gaussian.
347 @type alpha: float
348
349 @return: an array (length = 2*len(|inputSeries|) - 1) containing the smoothed signal.
350 @rtype: NumPy array
351 """
352
353
354
355
356 outputSeries = N.zeros((2*len(inputSeries) - 2,), N.Float)
357
358
359 res = inputSeries*N.exp(-0.5*(alpha*N.arange(len(inputSeries))/(len(inputSeries)-1))**2)
360
361
362 outputSeries[:len(inputSeries)] = res
363 outputSeries[len(inputSeries):] = res[-2:0:-1]
364
365 return outputSeries
366
368 """Returns n!
369
370 @param n: the n of n!.
371 @type: integer
372
373 @return: n!.
374 @rtype: integer
375 """
376 if n < 2:
377 return 1
378 else:
379
380 return N.multiply.reduce(N.arange(2,n+1))
381
383 """Intermediate function used to setup the calculation of spherical harmonics."""
384
385
386 c2 = N.sqrt((factorial(j+m) * factorial(j-m))) * factorial(j)
387
388 c3 = j+m
389 pp = []
390 aa = []
391 for p in range(0,j+1,1):
392 if (p <= c3) and (p >= m):
393
394 p1 = c3-p
395
396 p2 = j-p
397 p3 = p-m
398
399
400 a1 = ((-1)**p)*c2/(factorial(p1)*factorial(p2)*factorial(p)*factorial(p3))
401
402 pp.append((p1,p2,p3,p))
403
404 aa.append(a1)
405
406 return N.array(pp), N.array(aa)
407
409 """This function returns the r, theta and phi spherical coordinates from the x, y z cartesian coordinates.
410
411 @param x: the cartesian x.
412 @type x: float
413
414 @param y: the cartesian y.
415 @type y: float
416
417 @param z: the cartesian z.
418 @type z: float
419
420 @return: the r, theta and phi spherical coordinates..
421 @rtype: a list of three floats
422 """
423
424
425 r = N.sqrt(x**2 + y**2 + z**2)
426
427
428 theta = N.arccos(z/r)
429
430
431 phi = N.arctan2(y,x)
432
433 return r, theta, phi
434
436 """This function return the coordinates
437
438 @param pt: the coordinates of the point in the old basis.
439 @type pt: Scientific Vector
440
441 @param op: the coordinates of the new origin in the old basis.
442 @type op: Scientific Vector
443
444 @param ip: the coordinates of the new x axis in the old basis.
445 @type ip: Scientific Vector
446
447 @param jp: the coordinates of the new y axis in the old basis.
448 @type jp: Scientific Vector
449
450 @param kp: the coordinates of the new z axis in the old basis.
451 @type kp: Scientific Vector
452
453 @return: the coordinates of the point in the new basis.
454 @rtype: Scientific Vector
455 """
456
457 mat = Tensor(N.array([ip,jp,kp]).transpose())
458 matinv = mat.inverse()
459
460 mprim = matinv * (pt - op)
461
462 return mprim
463
464
465
466
467
469 """Returns the basis vectors for the simulation cell from the six crystallographic parameters.
470
471 @param parameters: a list of six floats defining the simulation cell geometry.
472 @type: parameters: list
473
474 @return: a list of three Scientific.Geometry.Vector objects representing respectively a, b and c
475 basis vectors.
476 @rtype: list
477 """
478
479 a, b, c, alpha, beta, gamma = parameters
480 e1 = Vector(a, 0.0, 0.0)
481 e2 = b*Vector(N.cos(gamma), N.sin(gamma), 0.0)
482 e3_x = N.cos(beta)
483 e3_y = (N.cos(alpha) - N.cos(beta)*N.cos(gamma)) / N.sin(gamma)
484 e3_z = N.sqrt(1.0 - e3_x**2 - e3_y**2)
485 e3 = c*Vector(e3_x, e3_y, e3_z)
486
487 return (e1, e2, e3)
488
490 """Returns a vector drawn from an uniform distribution within a circle of radius |r| and
491 orthogonal to vector |axis|.
492
493 @param r: float specifying the radius of the circle.
494 @type r: float
495
496 @param axis: the axis orthogonal to the plane where the vectors have to be generated.
497 @type axis: Scientific.Geometry.Vector object
498
499 @return: a vector pointing to a random point of the circle.
500 @rtype: Scientific.Geometry.Vector object
501 """
502 rsq = r*r
503 while 1:
504 x = N.array([uniform(-r, r), uniform(-r,r), uniform(-r,r)])
505 y = Vector(x) - (Vector(x)*axis/axis.length()**2)*axis
506 if y*y < rsq :
507 break
508 return y
509
511 """Returns a normalized vector drawn from an uniform distribution on the surface of a unit circle on a
512 plane orthogonal to |axis|.
513
514 @param axis: the axis orthogonal to the plane where the vectors have to be generated.
515 @type axis: Scientific.Geometry.Vector object
516
517 @return: A normalized vector defined in a unit disk orthogonal to |axis|
518 @rtype: Scientific.Geometry.Vector object
519 """
520 r = randomPointInCircle(1.0,axis)
521 return r.normal()
522
524 """Returns a normalized random vector on a plane or in space.
525
526 @param directions: if not None, a list of 2 Scientific.Vector that will define the plane on which
527 the vector should be generated.
528 @type directions: list of 2 Scientific.Vector or None
529
530 @return: a normalized random vector on a plane defined by |directions| or in space (|directions| = None).
531 @rtype: Scientific.Geometry.Vector object
532 """
533
534 if directions is None:
535 return randomDirection()
536
537 if len(directions) == 1:
538 return directions[0]
539
540 elif len(directions) == 2:
541
542 v1, v2 = directions
543
544
545 axis = (v1.cross(v2))
546
547 if axis.length() == 0.0:
548 raise Error('Your basis vectors are colinear. Can not build a plane out of them.')
549
550 return randomDirection2D(axis)
551