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