tools/i-pi/ipi/utils/nmtransform.py

   1 """Contains functions for doing the inverse and forward normal mode transforms.
   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 Classes:
  20    nm_trans: Uses matrix multiplication to do normal mode transformations.
  21    nm_rescale: Uses matrix multiplication to do ring polymer contraction
  22       or expansion.
  23    nm_fft: Uses fast-Fourier transforms to do normal modes transformations.
  24
  25 Functions:
  26    mk_nm_matrix: Makes a matrix to transform between the normal mode and bead
  27       representations.
  28    mk_rs_matrix: Makes a matrix to transform between one number of beads and
  29       another. Higher normal modes in the case of an expansion are set to zero.
  30 """
  31
  32 __all__ = ['nm_trans', 'nm_rescale', 'nm_fft']
  33
  34 import numpy as np
  35 from ipi.utils.messages import verbosity, info
  36
  37 def mk_nm_matrix(nbeads):
  38    """Gets the matrix that transforms from the bead representation
  39    to the normal mode representation.
  40
  41    If we return from this function a matrix C, then we transform between the
  42    bead and normal mode representation using q_nm = C . q_b, q_b = C.T . q_nm
  43
  44    Args:
  45       nbeads: The number of beads.
  46    """
  47
  48    b2nm = np.zeros((nbeads,nbeads))
  49    b2nm[0,:] = np.sqrt(1.0)
  50    for j in range(nbeads):
  51       for i in range(1, nbeads/2+1):
  52          b2nm[i,j] = np.sqrt(2.0)*np.cos(2*np.pi*j*i/float(nbeads))
  53       for i in range(nbeads/2+1, nbeads):
  54          b2nm[i,j] = np.sqrt(2.0)*np.sin(2*np.pi*j*i/float(nbeads))
  55    if (nbeads%2) == 0:
  56       b2nm[nbeads/2,0:nbeads:2] = 1.0
  57       b2nm[nbeads/2,1:nbeads:2] = -1.0
  58    return b2nm/np.sqrt(nbeads)
  59
  60 def mk_rs_matrix(nb1, nb2):
  61    """Gets the matrix that transforms a path with nb1 beads into one with
  62    nb2 beads.
  63
  64    If we return from this function a matrix T, then we transform between the
  65    system with nb1 bead and the system of nb2 beads using q_2 = T . q_1
  66
  67    Args:
  68       nb1: The initial number of beads.
  69       nb2: The final number of beads.
  70    """
  71
  72    if (nb1 == nb2):
  73       return np.identity(nb1,float)
  74    elif (nb1 > nb2):
  75       b1_nm = mk_nm_matrix(nb1)
  76       nm_b2 = mk_nm_matrix(nb2).T
  77
  78       #builds the "reduction" matrix that picks the normal modes we want to keep
  79       b1_b2 = np.zeros((nb2, nb1), float)
  80       b1_b2[0,0] = 1.0
  81       for i in range(1, nb2/2+1):
  82          b1_b2[i,i] = 1.0
  83          b1_b2[nb2-i, nb1-i] = 1.0
  84       if (nb2 % 2 == 0):
  85          #if we are contracting down to an even number of beads, then we have to
  86          #pick just one of the last degenerate modes to match onto the single
  87          #stiffest mode in the new path
  88          b1_b2[nb2/2, nb1-nb2/2] = 0.0
  89
  90       rs_b1_b2 = np.dot(nm_b2, np.dot(b1_b2, b1_nm))
  91       return rs_b1_b2*np.sqrt(float(nb2)/float(nb1))
  92    else:
  93       return mk_rs_matrix(nb2, nb1).T*(float(nb2)/float(nb1))
  94
  95
  96 class nm_trans:
  97    """Helper class to perform beads <--> normal modes transformation.
  98
  99    Attributes:
 100       _b2nm: The matrix to transform between the bead and normal mode
 101          representations.
 102       _nm2b: The matrix to transform between the normal mode and bead
 103          representations.
 104    """
 105
 106    def __init__(self, nbeads):
 107       """Initializes nm_trans.
 108
 109       Args:
 110          nbeads: The number of beads.
 111       """
 112
 113       self._b2nm = mk_nm_matrix(nbeads)
 114       self._nm2b = self._b2nm.T
 115
 116    def b2nm(self, q):
 117       """Transforms a matrix to the normal mode representation.
 118
 119       Args:
 120          q: A matrix with nbeads rows, in the bead representation.
 121       """
 122
 123       return np.dot(self._b2nm,q)
 124
 125    def nm2b(self, q):
 126       """Transforms a matrix to the bead representation.
 127
 128       Args:
 129          q: A matrix with nbeads rows, in the normal mode representation.
 130       """
 131
 132       return np.dot(self._nm2b,q)
 133
 134
 135 class nm_rescale:
 136    """Helper class to rescale a ring polymer between different number of beads.
 137
 138    Attributes:
 139       _b1tob2: The matrix to transform between a ring polymer with 'nbeads1'
 140          beads and another with 'nbeads2' beads.
 141       _b2tob1: The matrix to transform between a ring polymer with 'nbeads2'
 142          beads and another with 'nbeads1' beads.
 143    """
 144
 145    def __init__(self, nbeads1, nbeads2):
 146       """Initializes nm_rescale.
 147
 148       Args:
 149          nbeads1: The initial number of beads.
 150          nbeads2: The rescaled number of beads.
 151       """
 152
 153       self._b1tob2 = mk_rs_matrix(nbeads1,nbeads2)
 154       self._b2tob1 = self._b1tob2.T*(float(nbeads1)/float(nbeads2))
 155
 156    def b1tob2(self, q):
 157       """Transforms a matrix from one value of beads to another.
 158
 159       Args:
 160          q: A matrix with nbeads1 rows, in the bead representation.
 161       """
 162
 163       return np.dot(self._b1tob2,q)
 164
 165    def b2tob1(self, q):
 166       """Transforms a matrix from one value of beads to another.
 167
 168       Args:
 169          q: A matrix with nbeads2 rows, in the bead representation.
 170       """
 171
 172       return np.dot(self._b2tob1,q)
 173
 174
 175
 176 class nm_fft:
 177    """Helper class to perform beads <--> normal modes transformation
 178       using Fast Fourier transforms.
 179
 180    Attributes:
 181       fft: The fast-Fourier transform function to transform between the
 182          bead and normal mode representations.
 183       ifft: The inverse fast-Fourier transform function to transform
 184          between the normal mode and bead representations.
 185       qdummy: A matrix to hold a copy of the bead positions to transform
 186          them to the normal mode representation.
 187       qnmdummy: A matrix to hold a copy of the normal modes to transform
 188          them to the bead representation.
 189       nbeads: The number of beads.
 190       natoms: The number of atoms.
 191    """
 192
 193    def __init__(self, nbeads, natoms):
 194       """Initializes nm_trans.
 195
 196       Args:
 197          nbeads: The number of beads.
 198          natoms: The number of atoms.
 199       """
 200
 201       self.nbeads = nbeads
 202       self.natoms = natoms
 203       try:
 204          import pyfftw
 205          info("Import of PyFFTW successful", verbosity.medium)
 206          self.qdummy = pyfftw.n_byte_align_empty((nbeads, 3*natoms), 16, 'float32')
 207          self.qnmdummy = pyfftw.n_byte_align_empty((nbeads//2+1, 3*natoms), 16, 'complex64')
 208          self.fft = pyfftw.FFTW(self.qdummy, self.qnmdummy, axes=(0,), direction='FFTW_FORWARD')
 209          self.ifft = pyfftw.FFTW(self.qnmdummy, self.qdummy, axes=(0,), direction='FFTW_BACKWARD')
 210       except ImportError: #Uses standard numpy fft library if nothing better
 211                           #is available
 212          info("Import of PyFFTW unsuccessful, using NumPy library instead", verbosity.medium)
 213          self.qdummy = np.zeros((nbeads,3*natoms), dtype='float32')
 214          self.qnmdummy = np.zeros((nbeads//2+1,3*natoms), dtype='complex64')
 215          def dummy_fft(self):
 216             self.qnmdummy = np.fft.rfft(self.qdummy, axis=0)
 217          def dummy_ifft(self):
 218             self.qdummy = np.fft.irfft(self.qnmdummy, n=self.nbeads, axis=0)
 219          self.fft = lambda: dummy_fft(self)
 220          self.ifft = lambda: dummy_ifft(self)
 221
 222    def b2nm(self, q):
 223       """Transforms a matrix to the normal mode representation.
 224
 225       Args:
 226          q: A matrix with nbeads rows and 3*natoms columns,
 227             in the bead representation.
 228       """
 229
 230       if self.nbeads == 1:
 231          return q
 232       self.qdummy[:] = q
 233       self.fft()
 234       if self.nbeads == 2:
 235          return self.qnmdummy.real/np.sqrt(self.nbeads)
 236
 237       nmodes = self.nbeads/2
 238
 239       self.qnmdummy /= np.sqrt(self.nbeads)
 240       qnm = np.zeros(q.shape)
 241       qnm[0,:] = self.qnmdummy[0,:].real
 242
 243       if self.nbeads % 2 == 0:
 244          self.qnmdummy[1:-1,:] *= np.sqrt(2)
 245          (qnm[1:nmodes,:], qnm[self.nbeads:nmodes:-1,:]) = (self.qnmdummy[1:-1,:].real, self.qnmdummy[1:-1,:].imag)
 246          qnm[nmodes,:] = self.qnmdummy[nmodes,:].real
 247       else:
 248          self.qnmdummy[1:,:] *= np.sqrt(2)
 249          (qnm[1:nmodes+1,:], qnm[self.nbeads:nmodes:-1,:]) = (self.qnmdummy[1:,:].real, self.qnmdummy[1:,:].imag)
 250
 251       return qnm
 252
 253    def nm2b(self, qnm):
 254       """Transforms a matrix to the bead representation.
 255
 256       Args:
 257          qnm: A matrix with nbeads rows and 3*natoms columns,
 258             in the normal mode representation.
 259       """
 260
 261       if self.nbeads == 1:
 262          return qnm
 263       if self.nbeads == 2:
 264          self.qnmdummy[:] = qnm
 265          self.ifft()
 266          return self.qdummy*np.sqrt(self.nbeads)
 267
 268       nmodes = self.nbeads/2
 269       odd = self.nbeads - 2*nmodes  # 0 if even, 1 if odd
 270
 271       qnm_complex = np.zeros((nmodes+1, len(qnm[0,:])), complex)
 272       qnm_complex[0,:] = qnm[0,:]
 273       if not odd:
 274          (qnm_complex[1:-1,:].real, qnm_complex[1:-1,:].imag) = (qnm[1:nmodes,:], qnm[self.nbeads:nmodes:-1,:])
 275          qnm_complex[1:-1,:] /= np.sqrt(2)
 276          qnm_complex[nmodes,:] = qnm[nmodes,:]
 277       else:
 278          (qnm_complex[1:,:].real, qnm_complex[1:,:].imag) = (qnm[1:nmodes+1,:], qnm[self.nbeads:nmodes:-1,:])
 279          qnm_complex[1:,:] /= np.sqrt(2)
 280
 281       self.qnmdummy[:] = qnm_complex
 282       self.ifft()
 283       return self.qdummy*np.sqrt(self.nbeads)