← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/i-pi/ipi/utils/mathtools.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/i-pi/ipi/utils/mathtools.py
 
1
"""Contains simple algorithms.
 
2
 
 
3
Copyright (C) 2013, Joshua More and Michele Ceriotti
 
4
 
 
5
This program is free software: you can redistribute it and/or modify
 
6
it under the terms of the GNU General Public License as published by
 
7
the Free Software Foundation, either version 3 of the License, or
 
8
(at your option) any later version.
 
9
 
 
10
This program is distributed in the hope that it will be useful,
 
11
but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 
13
GNU General Public License for more details.
 
14
 
 
15
You should have received a copy of the GNU General Public License
 
16
along with this program. If not, see <http.//www.gnu.org/licenses/>.
 
17
 
 
18
 
 
19
Functions:
 
20
   matrix_exp: Computes the exponential of a square matrix via a Taylor series.
 
21
   stab_cholesky: A numerically stable version of the Cholesky decomposition.
 
22
   h2abc: Takes the representation of the system box in terms of an upper
 
23
      triangular matrix of column vectors, and returns the representation in
 
24
      terms of the lattice vector lengths and the angles between them
 
25
      in radians.
 
26
   h2abc_deg: Takes the representation of the system box in terms of an upper
 
27
      triangular matrix of column vectors, and returns the representation in
 
28
      terms of the lattice vector lengths and the angles between them in
 
29
      degrees.
 
30
   abc2h: Takes the representation of the system box in terms of the lattice
 
31
      vector lengths and the angles between them, and returns the
 
32
      representation in terms of an upper triangular lattice vector matrix.
 
33
   invert_ut3x3: Inverts a 3*3 upper triangular matrix.
 
34
   det_ut3x3(h): Finds the determinant of a 3*3 upper triangular matrix.
 
35
   eigensystem_ut3x3: Finds the eigenvector matrix and eigenvalues of a 3*3
 
36
      upper triangular matrix
 
37
   exp_ut3x3: Computes the exponential of a 3*3 upper triangular matrix.
 
38
   root_herm: Computes the square root of a positive-definite hermitian
 
39
      matrix.
 
40
   logsumlog: Routine to accumulate the logarithm of a sum
 
41
"""
 
42
 
 
43
__all__ = ['matrix_exp', 'stab_cholesky', 'h2abc', 'h2abc_deg', 'abc2h',
 
44
           'invert_ut3x3', 'det_ut3x3', 'eigensystem_ut3x3', 'exp_ut3x3',
 
45
            'root_herm', 'logsumlog' ]
 
46
 
 
47
import numpy as np
 
48
import math
 
49
from ipi.utils.messages import verbosity, warning
 
50
 
 
51
def logsumlog(lasa, lbsb):
 
52
   """Computes log(|A+B|) and sign(A+B) given log(|A|), log(|B|),
 
53
   sign(A), sign(B).
 
54
 
 
55
   Args:
 
56
      lasa: (log(|A|), sign(A)) as a tuple
 
57
      lbsb: (log(|B|), sign(B)) as a tuple
 
58
 
 
59
   Returns:
 
60
      (log(|A+B|), sign(A+B)) as a tuple
 
61
   """
 
62
 
 
63
   (la,sa) = lasa
 
64
   (lb,sb) = lbsb
 
65
 
 
66
   if (la > lb):
 
67
      sr = sa
 
68
      lr = la + np.log(1.0 + sb*np.exp(lb-la))
 
69
   else:
 
70
      sr = sb
 
71
      lr = lb + np.log(1.0 + sa*np.exp(la-lb))
 
72
 
 
73
   return (lr,sr)
 
74
 
 
75
def matrix_exp(M, ntaylor=15, nsquare=15):
 
76
   """Computes the exponential of a square matrix via a Taylor series.
 
77
 
 
78
   Calculates the matrix exponential by first calculating exp(M/(2**nsquare)),
 
79
   then squaring the result the appropriate number of times.
 
80
 
 
81
   Args:
 
82
      M: Matrix to be exponentiated.
 
83
      ntaylor: Optional integer giving the number of terms in the Taylor series.
 
84
         Defaults to 15.
 
85
      nsquare: Optional integer giving how many times the original matrix will
 
86
         be halved. Defaults to 15.
 
87
 
 
88
   Returns:
 
89
      The matrix exponential of M.
 
90
   """
 
91
 
 
92
   n = M.shape[1]
 
93
   tc = np.zeros(ntaylor+1)
 
94
   tc[0] = 1.0
 
95
   for i in range(ntaylor):
 
96
      tc[i+1] = tc[i]/(i+1)
 
97
 
 
98
   SM = np.copy(M)/2.0**nsquare
 
99
 
 
100
   EM = np.identity(n,float)*tc[ntaylor]
 
101
   for i in range(ntaylor-1,-1,-1):
 
102
      EM = np.dot(SM,EM)
 
103
      EM += np.identity(n)*tc[i]
 
104
 
 
105
   for i in range(nsquare):
 
106
      EM = np.dot(EM,EM)
 
107
   return EM
 
108
 
 
109
def stab_cholesky(M):
 
110
   """ A numerically stable version of the Cholesky decomposition.
 
111
 
 
112
   Used in the GLE implementation. Since many of the matrices used in this
 
113
   algorithm have very large and very small numbers in at once, to handle a
 
114
   wide range of frequencies, a naive algorithm can end up having to calculate
 
115
   the square root of a negative number, which breaks the algorithm. This is
 
116
   due to numerical precision errors turning a very tiny positive eigenvalue
 
117
   into a tiny negative value.
 
118
 
 
119
   Instead of this, an LDU decomposition is used, and any small negative numbers
 
120
   in the diagonal D matrix are assumed to be due to numerical precision errors,
 
121
   and so are replaced with zero.
 
122
 
 
123
   Args:
 
124
      M: The matrix to be decomposed.
 
125
   """
 
126
 
 
127
   n = M.shape[1]
 
128
   D = np.zeros(n,float)
 
129
   L = np.zeros(M.shape,float)
 
130
   for i in range(n):
 
131
      L[i,i] = 1.
 
132
      for j in range(i):
 
133
         L[i,j] = M[i,j]
 
134
         for k in range(j):
 
135
            L[i,j] -= L[i,k]*L[j,k]*D[k]
 
136
         if (not D[j] == 0.0):
 
137
            L[i,j] = L[i,j]/D[j]
 
138
      D[i] = M[i,i]
 
139
      for k in range(i):
 
140
         D[i] -= L[i,k]*L[i,k]*D[k]
 
141
 
 
142
   S = np.zeros(M.shape,float)
 
143
   for i in range(n):
 
144
      if (D[i]>0):
 
145
         D[i] = math.sqrt(D[i])
 
146
      else:
 
147
         warning("Zeroing negative element in stab-cholesky decomposition: " + str(D[i]), verbosity.low)
 
148
         D[i] = 0
 
149
      for j in range(i+1):
 
150
         S[i,j] += L[i,j]*D[j]
 
151
   return S
 
152
 
 
153
def h2abc(h):
 
154
   """Returns a description of the cell in terms of the length of the
 
155
      lattice vectors and the angles between them in radians.
 
156
 
 
157
   Args:
 
158
      h: Cell matrix in upper triangular column vector form.
 
159
 
 
160
   Returns:
 
161
      A list containing the lattice vector lengths and the angles between them.
 
162
   """
 
163
 
 
164
   a = float(h[0,0])
 
165
   b = math.sqrt(h[0,1]**2 + h[1,1]**2)
 
166
   c = math.sqrt(h[0,2]**2 + h[1,2]**2 + h[2,2]**2)
 
167
   gamma = math.acos(h[0,1]/b)
 
168
   beta = math.acos(h[0,2]/c)
 
169
   alpha = math.acos(np.dot(h[:,1], h[:,2])/(b*c))
 
170
 
 
171
   return a, b, c, alpha, beta, gamma
 
172
 
 
173
def h2abc_deg(h):
 
174
   """Returns a description of the cell in terms of the length of the
 
175
      lattice vectors and the angles between them in degrees.
 
176
 
 
177
   Args:
 
178
      h: Cell matrix in upper triangular column vector form.
 
179
 
 
180
   Returns:
 
181
      A list containing the lattice vector lengths and the angles between them
 
182
      in degrees.
 
183
   """
 
184
 
 
185
   (a, b, c, alpha, beta, gamma) = h2abc(h)
 
186
   return a, b, c, alpha*180/math.pi, beta*180/math.pi, gamma*180/math.pi
 
187
 
 
188
def abc2h(a, b, c, alpha, beta, gamma):
 
189
   """Returns a lattice vector matrix given a description in terms of the
 
190
   lattice vector lengths and the angles in between.
 
191
 
 
192
   Args:
 
193
      a: First cell vector length.
 
194
      b: Second cell vector length.
 
195
      c: Third cell vector length.
 
196
      alpha: Angle between sides b and c in radians.
 
197
      beta: Angle between sides a and c in radians.
 
198
      gamma: Angle between sides b and a in radians.
 
199
 
 
200
   Returns:
 
201
      An array giving the lattice vector matrix in upper triangular form.
 
202
   """
 
203
 
 
204
   h = np.zeros((3,3) ,float)
 
205
   h[0,0] = a
 
206
   h[0,1] = b*math.cos(gamma)
 
207
   h[0,2] = c*math.cos(beta)
 
208
   h[1,1] = b*math.sin(gamma)
 
209
   h[1,2] = (b*c*math.cos(alpha) - h[0,1]*h[0,2])/h[1,1]
 
210
   h[2,2] = math.sqrt(c**2 - h[0,2]**2 - h[1,2]**2)
 
211
   return h
 
212
 
 
213
def invert_ut3x3(h):
 
214
   """Inverts a 3*3 upper triangular matrix.
 
215
 
 
216
   Args:
 
217
      h: An upper triangular 3*3 matrix.
 
218
 
 
219
   Returns:
 
220
      The inverse matrix of h.
 
221
   """
 
222
 
 
223
   ih = np.zeros((3,3), float)
 
224
   for i in range(3):
 
225
      ih[i,i] = 1.0/h[i,i]
 
226
   ih[0,1] = -ih[0,0]*h[0,1]*ih[1,1]
 
227
   ih[1,2] = -ih[1,1]*h[1,2]*ih[2,2]
 
228
   ih[0,2] = -ih[1,2]*h[0,1]*ih[0,0] - ih[0,0]*h[0,2]*ih[2,2]
 
229
   return ih
 
230
 
 
231
def eigensystem_ut3x3(p):
 
232
   """Finds the eigenvector matrix of a 3*3 upper-triangular matrix.
 
233
 
 
234
   Args:
 
235
      p: An upper triangular 3*3 matrix.
 
236
 
 
237
   Returns:
 
238
      An array giving the 3 eigenvalues of p, and the eigenvector matrix of p.
 
239
   """
 
240
 
 
241
   eigp = np.zeros((3,3), float)
 
242
   eigvals = np.zeros(3, float)
 
243
 
 
244
   for i in range(3):
 
245
      eigp[i,i] = 1
 
246
   eigp[0,1] = -p[0,1]/(p[0,0] - p[1,1])
 
247
   eigp[1,2] = -p[1,2]/(p[1,1] - p[2,2])
 
248
   eigp[0,2] = -(p[0,1]*p[1,2] - p[0,2]*p[1,1] + p[0,2]*p[2,2])/((p[0,0] - p[2,2])*(p[2,2] - p[1,1]))
 
249
 
 
250
   for i in range(3):
 
251
      eigvals[i] = p[i,i]
 
252
   return eigp, eigvals
 
253
 
 
254
def det_ut3x3(h):
 
255
   """Calculates the determinant of a 3*3 upper triangular matrix.
 
256
 
 
257
   Note that the volume of the system box when the lattice vector matrix is
 
258
   expressed as a 3*3 upper triangular matrix is given by the determinant of
 
259
   this matrix.
 
260
 
 
261
   Args:
 
262
      h: An upper triangular 3*3 matrix.
 
263
 
 
264
   Returns:
 
265
      The determinant of h.
 
266
   """
 
267
 
 
268
   return h[0,0]*h[1,1]*h[2,2]
 
269
 
 
270
MINSERIES=1e-8
 
271
def exp_ut3x3(h):
 
272
   """Computes the matrix exponential for a 3x3 upper-triangular matrix.
 
273
 
 
274
   Note that for 3*3 upper triangular matrices this is the best method, as
 
275
   it is stable. This is terms which become unstable as the
 
276
   denominator tends to zero are calculated via a Taylor series in this limit.
 
277
 
 
278
   Args:
 
279
      h: An upper triangular 3*3 matrix.
 
280
 
 
281
   Returns:
 
282
      The matrix exponential of h.
 
283
   """
 
284
   eh = np.zeros((3,3), float)
 
285
   e00 = math.exp(h[0,0])
 
286
   e11 = math.exp(h[1,1])
 
287
   e22 = math.exp(h[2,2])
 
288
   eh[0,0] = e00
 
289
   eh[1,1] = e11
 
290
   eh[2,2] = e22
 
291
 
 
292
   if (abs((h[0,0] - h[1,1])/h[0,0])>MINSERIES):
 
293
      r01 = (e00 - e11)/(h[0,0] - h[1,1])
 
294
   else:
 
295
      r01 = e00*(1 + (h[0,0] - h[1,1])*(0.5 + (h[0,0] - h[1,1])/6.0))
 
296
   if (abs((h[1,1] - h[2,2])/h[1,1])>MINSERIES):
 
297
      r12 = (e11 - e22)/(h[1,1] - h[2,2])
 
298
   else:
 
299
      r12 = e11*(1 + (h[1,1] - h[2,2])*(0.5 + (h[1,1] - h[2,2])/6.0))
 
300
   if (abs((h[2,2] - h[0,0])/h[2,2])>MINSERIES):
 
301
      r02 = (e22 - e00)/(h[2,2] - h[0,0])
 
302
   else:
 
303
      r02 = e22*(1 + (h[2,2] - h[0,0])*(0.5 + (h[2,2] - h[0,0])/6.0))
 
304
 
 
305
   eh[0,1] = h[0,1]*r01
 
306
   eh[1,2] = h[1,2]*r12
 
307
 
 
308
   eh[0,2] = h[0,2]*r02
 
309
   if (abs((h[2,2] - h[0,0])/h[2,2])>MINSERIES):
 
310
      eh[0,2] += h[0,1]*h[0,2]*(r01 - r12)/(h[0,0] - h[2,2])
 
311
   elif (abs((h[1,1] - h[0,0])/h[1,1])>MINSERIES):
 
312
      eh[0,2] += h[0,1]*h[0,2]*(r12 - r02)/(h[1,1] - h[0,0])
 
313
   elif (abs((h[1,1]-h[2,2])/h[1,1])>MINSERIES):
 
314
      eh[0,2] += h[0,1]*h[0,2]*(r02 - r01)/(h[2,2] - h[1,1])
 
315
   else:
 
316
      eh[0,2] += h[0,1]*h[0,2]*e00/24.0*(12.0 + 4*(h[1,1] + h[2,2] - 2*h[0,0]) + (h[1,1] - h[0,0])*(h[2,2] - h[0,0]))
 
317
 
 
318
   return eh
 
319
 
 
320
def root_herm(A):
 
321
   """Gives the square root of a hermitian matrix with real eigenvalues.
 
322
 
 
323
   Args:
 
324
      A: A Hermitian matrix.
 
325
 
 
326
   Returns:
 
327
      A matrix such that itself matrix multiplied by its transpose gives the
 
328
      original matrix.
 
329
   """
 
330
 
 
331
   if not (abs(A.T - A) < 1e-10).all():
 
332
      raise ValueError("Non-Hermitian matrix passed to root_herm function")
 
333
   eigvals, eigvecs = np.linalg.eigh(A)
 
334
   ndgrs = len(eigvals)
 
335
   diag = np.zeros((ndgrs,ndgrs))
 
336
   for i in range(ndgrs):
 
337
      if eigvals[i] >= 0:
 
338
         diag[i,i] = math.sqrt(eigvals[i])
 
339
      else:
 
340
         warning("Zeroing negative element in matrix square root: " + str(eigvals[i]), verbosity.low)
 
341
         diag[i,i] = 0
 
342
   return np.dot(eigvecs, np.dot(diag, eigvecs.T))
 
343