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