1"""Contains simple algorithms.
2
3Copyright (C) 2013, Joshua More and Michele Ceriotti
4
5This program is free software: you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation, either version 3 of the License, or
8(at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with this program. If not, see <http.//www.gnu.org/licenses/>.
17
18
19Functions:
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
47import numpy as np
48import math
49from ipi.utils.messages import verbosity, warning
50
51def 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
75def 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
109def 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
153def 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
173def 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
188def 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
213def 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
231def 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
254def 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
270MINSERIES=1e-8
271def 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
320def 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
344