Source code for Inelastica.physics.mesh

"""

:mod:`Inelastica.physics.mesh`
==============================

Written by Thomas Frederiksen.

Classes
-------

.. autosummary::
   :toctree:

   kmesh

.. currentmodule:: Inelastica.physics.mesh

"""
from __future__ import print_function

import Inelastica.math as MM
import numpy as N

# This function is used in the class below for each of the
# three spatial directions (in momentum space)


[docs] def generatelinmesh(Nk): "Generate Nk points sampling linearly the interval [-0.5;0.5]" kpts = [(ii*1.0+0.5)/Nk-0.5 for ii in range(Nk)] wgts = [1.0/Nk for ii in range(Nk)] return N.array(kpts), N.array(wgts)
[docs] class kmesh(object): """ Create a k-mesh samling where each of the three components sample the range [-0.5,0.5]). They are not in reciprocal space. Instead they correspond to the mathematical orthogonal space that is Fourier transformed. Variables: k : Array of 3D k-points (shape = [NNk,3]) w : Array of integration weights (shape = [4,NNk]), where w[0] are the full weights and w[i=1,2,3] are weights with a reduced number of GK points for the three axes. These are used for error estimates within the GK method. NNk : Total number of k-points explicitly in the mesh Nk : Array with sampling for each axis type : List of sampling type for each axis """
[docs] def __init__(self, Nk1, Nk2, Nk3, meshtype=['LIN', 'LIN', 'LIN'], invsymmetry=True): """ Returns an instance of a k-mesh with each of the k-vector axes sampled either linearly (LIN) or using a Gauss-Kronrod (GK) scheme. An axis i sampled by (Nki,LIN) generates Nk points, while (Nki,GK) returns 2*Nk+1 points. By applying inversion symmetry (k=-k) the number of points in the mesh can be reduced by a factor two (except when the Gamma point is included which have no partner). """ self.Nk = N.array([Nk1, Nk2, Nk3]) self.type = meshtype self.genkmesh() self.invsymmetry = invsymmetry if invsymmetry: self.SymmetryReduce()
[docs] def genkmesh(self): "Generate mesh for a given sampling of the axes" self.k = [] self.w = [] errorw = [] for i in range(3): # loop over the three k-components if self.type[i].upper() == 'GK' or self.type[i].upper() == 'GAUSSKRONROD': self.type[i] = 'GK' if self.Nk[i] > 1: # GK-method fails with fewer points kpts, wgts, ew = MM.GaussKronrod(self.Nk[i]) self.k.append(kpts) self.w.append(wgts) errorw.append(ew) else: print('Kmesh.py: GK method requires Nk=%i>1'%(self.Nk[i])) sys.exit(1) elif self.type[i].upper() == 'LIN' or self.type[i].upper() == 'LINEAR': self.type[i] = 'LIN' kpts, wgts = generatelinmesh(self.Nk[i]) self.k.append(kpts) self.w.append(wgts) errorw.append(wgts) else: print('Kmesh.py: Unknown meshtype:', self.type[i].upper()) self.Nk[i] = len(self.k[i]) self.NNk = N.prod(self.Nk) print('Kmesh.py: Generating mesh:') print(' ... type = ', self.type) print(' ... Nk = ', self.Nk) # repete out in 3D kpts = N.zeros((self.NNk, 3)) # Array of k-points wgts = N.ones((4, self.NNk)) # (wgts, errorw1, errorw2, errorw3) nn = 0 for i in range(self.Nk[0]): for j in range(self.Nk[1]): for k in range(self.Nk[2]): kpts[nn, :] = [self.k[0][i], self.k[1][j], self.k[2][k]] wgts[0, nn] = self.w[0][i]*self.w[1][j]*self.w[2][k] wgts[1, nn] = errorw[0][i]*self.w[1][j]*self.w[2][k] wgts[2, nn] = self.w[0][i]*errorw[1][j]*self.w[2][k] wgts[3, nn] = self.w[0][i]*self.w[1][j]*errorw[2][k] nn += 1 print(' ... NNk = %i, sum(wgts) = %.8f'%(self.NNk, N.sum(wgts[0]))) print(' ... sum(errorw) = (%.8f,%.8f,%.8f)'%tuple(N.sum(wgts[i+1]) for i in range(3))) self.k = kpts self.w = wgts
[docs] def SymmetryReduce(self): """ Remove duplicates for symmetry INVERSION SYMMETRY: If the Bloch function \psi(k) = exp(ikr)u(k), with crystal momentum k, is an eigenstate of the Schroedinger equation then also \psi^\dagger(k) = exp(-ikr)u^\dagger(k) with crystal momentum -k, is an eigenstate with same eigenvalue. Hence E(k) = E(-k). TIME REVERSAL SYMMETRY: t,\psi(r,t) --> -t,\psi^\dagger(r,-t). T(k) = T(-k). (Elastic) propagation from L to R is always identical to propagation from R to L. """ print(' ... Applying inversion symmetry (the simple way)') # No brute force (and therefore terribly slow) pairing here indx = [[ii, 2] for ii in range(self.NNk // 2, self.NNk)] # Keep the last half of the k-points with double weight k0 = self.k[self.NNk // 2] if N.dot(k0, k0) == 0: # Gamma in the kptlist indx[0] = [self.NNk // 2, 1] # lower weight to one indx, weight = N.array([ii[0] for ii in indx]), N.array([ii[1] for ii in indx]) kpts, wgts = self.k[indx], self.w[:, indx]*weight self.k = kpts self.NNk = len(kpts) self.w = wgts print(' ... NNk = %i, sum(wgts) = %.8f'%(self.NNk, N.sum(wgts[0]))) print(' ... sum(errorw) = (%.8f,%.8f,%.8f)'%tuple(N.sum(wgts[i+1]) for i in range(3)))
[docs] def mesh2file(self, fn): "Writes the k-mesh to file" f = open(fn, 'w') f.write('# k1 k2 k3 w0 w1 w2 w3\n') for i in range(self.NNk): s = '' for j in range(3): s += '%.8f '%self.k[i, j] for j in range(4): s += '%.8e '%self.w[j, i] f.write(s+'\n') f.close()
[docs] def test(): """ Test function """ mesh = kmesh(4, 3, 3, meshtype=['LIN', 'GK', 'LIN'], invsymmetry=False) keys = list(mesh.__dict__.keys()) print(keys) print(mesh.Nk) print(mesh.NNk) print(mesh.type) print(mesh.invsymmetry) #print 'ki wi' #for i in range(len(mesh.k)): # print mesh.k[i], mesh.w[0, i] mesh.mesh2file('mesh-test.dat') print('Integrate some simple functions over [-0.5,0.5]:') print(' f(x,y,z)=1 => \int f dxdydz =', N.sum(mesh.w[0])) for i, s in enumerate(['x', 'y', 'z']): f = mesh.k[:, i] print(' f(x,y,z)=%s => \int f dxdydz ='%s, N.sum(f*mesh.w[0])) f = mesh.k[:, 0]*mesh.k[:, 1]*mesh.k[:, 2] print(' f(x,y,z)=x*y*z => \int f dxdydz =', N.sum(f*mesh.w[0])) f = (mesh.k[:, 0]+1)*(mesh.k[:, 1]+1)*(mesh.k[:, 2]+1) print(' f(x,y,z)=(x+1)*(y+1)*(z+1) => \int f dxdydz =', N.sum(f*mesh.w[0])) for i, s in enumerate(['x', 'y', 'z']): f = N.cos(2*mesh.k[:, i]) print(' f(x,y,z)=cos(2%s) => \int f dxdydz ='%s, N.sum(f*mesh.w[0]), '[exact: sin(1) ~ 0.841470984807897]')
if __name__ == '__main__': test()