TBit (empty) → 0.4.2.0
raw patch · 24 files changed
+1875/−0 lines, 24 filesdep +basedep +containersdep +deepseqsetup-changed
Dependencies added: base, containers, deepseq, fgl, free, haskeline, hmatrix, integration, list-extras, mtl, numeric-tools, parallel, stream-fusion
Files
- LICENSE +12/−0
- Setup.hs +2/−0
- TBit.cabal +84/−0
- TBit/Electronic/Conductance.hs +37/−0
- TBit/Framework.hs +139/−0
- TBit/Hamiltonian/Builder/Decompactification.hs +74/−0
- TBit/Hamiltonian/Builder/Examples.hs +93/−0
- TBit/Hamiltonian/Builder/Matrification.hs +51/−0
- TBit/Hamiltonian/Builder/PrimitiveLattice.hs +105/−0
- TBit/Hamiltonian/Builder/Terms.hs +138/−0
- TBit/Hamiltonian/Eigenstates.hs +41/−0
- TBit/Magnetic/OrbitalMagnetization.hs +144/−0
- TBit/Numerical/Derivative.hs +25/−0
- TBit/Numerical/Integration.hs +55/−0
- TBit/Parameterization.hs +88/−0
- TBit/Plots.hs +108/−0
- TBit/Systems/HoneycombLattice.hs +64/−0
- TBit/Systems/KagomeLattice.hs +132/−0
- TBit/Systems/SquareLattice.hs +20/−0
- TBit/Toolbox.hs +56/−0
- TBit/Topological/Chern.hs +110/−0
- TBit/Topological/Curvature.hs +83/−0
- TBit/Topological/Spillage.hs +49/−0
- TBit/Types.hs +165/−0
+ LICENSE view
@@ -0,0 +1,12 @@+Copyright (c) 2014, Matthew W. Daniels+All rights reserved.++Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ TBit.cabal view
@@ -0,0 +1,84 @@+-- Initial TightBinding.cabal generated by cabal init. For further +-- documentation, see http://haskell.org/cabal/users-guide/++-- The name of the package.+name: TBit++-- The package version. See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented.+-- http://www.haskell.org/haskellwiki/Package_versioning_policy+-- PVP summary: +-+------- breaking API changes+-- | | +----- non-breaking API additions+-- | | | +--- code changes with no API change+version: 0.4.2.0++-- A short (one-line) description of the package.+synopsis: Utilities for condensed matter physics tight binding calculations.++-- A longer description of the package.+description: TBit provides tools for parameterizing and computing condensed matter physics quantities based on tight-binding models. It provides utitilies for computing Chern numbers and Berry curvatures of electronic band structure, generating gnuplot-readable dispersion plots, and calculating assorted quantities such as orbital magnetization and Nernst conductivity.++-- The license under which the package is released.+license: BSD3++-- The file containing the license text.+license-file: LICENSE++-- The package author(s).+author: Matthew Daniels++-- An email address to which users can send suggestions, bug reports, and +-- patches.+maintainer: danielsmw@gmail.com++-- A copyright notice.+-- copyright: ++category: Science++build-type: Simple++-- Constraint on the version of Cabal needed to build this package.+cabal-version: >=1.8+++library+ -- Modules exported by the library.+ exposed-modules: TBit.Types+ , TBit.Plots+ , TBit.Framework+ , TBit.Parameterization+ , TBit.Toolbox++ , TBit.Systems.KagomeLattice+ , TBit.Systems.HoneycombLattice+ , TBit.Systems.SquareLattice++ , TBit.Topological.Chern+ , TBit.Topological.Curvature+ , TBit.Topological.Spillage++ , TBit.Magnetic.OrbitalMagnetization+ , TBit.Electronic.Conductance++ , TBit.Numerical.Derivative+ , TBit.Numerical.Integration++ , TBit.Hamiltonian.Eigenstates+ , TBit.Hamiltonian.Builder.Decompactification+ , TBit.Hamiltonian.Builder.PrimitiveLattice+ , TBit.Hamiltonian.Builder.Matrification+ , TBit.Hamiltonian.Builder.Terms+ , TBit.Hamiltonian.Builder.Examples++ extensions: DoAndIfThenElse, FlexibleInstances, FlexibleContexts, TypeFamilies+ + -- Modules included in this library but not exported.+ -- other-modules: + + -- Other library packages from which modules are imported.+ build-depends: base ==4.6.*, containers ==0.5.*,+ hmatrix ==0.16.*, stream-fusion ==0.1.*,+ mtl ==2.2.*, parallel ==3.2.*, deepseq ==1.3.*,+ integration ==0.2.*, numeric-tools ==0.2.*, fgl ==5.5.*,+ list-extras ==0.4.*, free ==4.9.*, haskeline ==0.7.*
+ TBit/Electronic/Conductance.hs view
@@ -0,0 +1,37 @@+module TBit.Electronic.Conductance (nernstConductivity) where++import Data.Complex (realPart)+import Control.Monad (liftM)+import Numeric.LinearAlgebra.HMatrix (size, vector)++import TBit.Types+import TBit.Framework+import TBit.Parameterization+import TBit.Hamiltonian.Eigenstates+import TBit.Topological.Curvature++nernstConductivity :: Hamiltonian -> Parameterized Double+nernstConductivity h = bzIntegral (nernstIntegrand h)++nernstIntegrand :: Hamiltonian -> Wavevector -> Parameterized Double+nernstIntegrand h k = liftM sum . sequence+ $ [ nernstIntegrand' h b k+ | b <- [0..pred dim] ]+ where dim = fst $ size $ h $ vector [0,0]++nernstIntegrand' :: Hamiltonian -> BandIndex -> Wavevector -> Parameterized Double+nernstIntegrand' h n k = + do omega <- bandCurvature n h k+ energy <- eigenenergies h k+ beta <- getScalar "beta"+ mu <- getScalar "mu"+ let f = f' (realPart $ beta) (realPart $ mu)+ let s = negate+ $ f (energy!!n) * (log.f) (energy!!n)+ + (1 - f (energy!!n)) * log (1 - f (energy!!n) + eps)+ return $ if (energy!!n) < (realPart mu)+ then omega * s+ else 0.0++ where f' beta mu e = 1.0 / (1.0 + (exp $ beta * (e - mu)))+ eps = 1.0e-50
+ TBit/Framework.hs view
@@ -0,0 +1,139 @@+{-# LANGUAGE NoImplicitPrelude #-}+module TBit.Framework ( meshBZ+ , meshBigBZ+ , bzIntegral+ , bzIntegral'+ , bzIntegral''+ , kPath ) where++import TBit.Types+import TBit.Parameterization+import TBit.Numerical.Integration++import Data.Foldable hiding (foldl', concatMap)+import Data.List.Stream hiding (sum)+import Data.Traversable (traverse)+import qualified Data.Map as M++import Numeric.LinearAlgebra.HMatrix (norm_2, scale, (!))++import Control.Monad (liftM)+import Control.Monad.State (get)+import Prelude.Listless ++{-|+ Given a list of points in /k/-space, return a list of+ points that interpolates affine paths between them, in+ turn; a typical usage case might be++ > kPath [gammaPoint, kPoint, mPoint, gammaPoint]++ which is used in the 'Plots.bandPlot' function. The spacing+ between points on the interpolated path is determined by+ the 'Types.meshingData' parameter.+-}+kPath :: [Wavevector] -> Parameterized ([Wavevector])+kPath [] = return []+kPath (k:[]) = return [k]+kPath (k:k':ks) = do (Spacing eps) <- getMesh+ nextPath <- kPath (k':ks)+ return $ [ k + scale n (dk eps)+ | n <- takeWhile ((>) (norm_2 (k'-k)) . (*) eps)+ [0.0, 1.0 ..]]+ ++ nextPath++ where dk eps = scale eps $ normalize (k'-k)+ normalize v = scale (1.0/(norm_2 v)) v++{-|+ Given a function defined on the Brillouin zone, evaluate it everywhere+ by doing nested single integrations and using Takahashi and Mori's+ Tanh-Sinh quadrature method. Should be robust against singularities+ and the like, and is /properly/ set up for massive parallelization (via+ EdwardKmett's integration library, certainly not mine). Watch Dirac run+ this thing. Only implemented in 2D.+-}+bzIntegral'' :: (Wavevector -> Parameterized Double) -> Parameterized Double+bzIntegral'' f = do (b1:b2:_) <- recipPrimitiveLattice+ ps <- get+ let g s t = f (scale s b1 + scale t b2) :: Parameterized Double+ return . (*) (jac b1 b2)+ . fst + $ integrate (\x y -> either e id + . flip crunch ps + $ g x y) (0.0,1.0) (0.0,1.0)+ where e err = error (show err)+ jac a b = abs $ (a!0) * (b!1) - (a!1) * (b!0)+++{-|+ As 'bzIntegral', but also gives the absolute error as the second value+ in the returned tuple.+-}+bzIntegral' :: (Wavevector -> Parameterized Double) -> Parameterized (Double, Double)+bzIntegral' f = do (b1:b2:_) <- recipPrimitiveLattice+ ps <- get+ let g s t = f (scale s b1 + scale t b2) :: Parameterized Double+ return . (\(x,y) -> (x * (jac b1 b2), y))+ $ integrate (\x y -> either e id + . flip crunch ps + $ g x y) (0.0,1.0) (0.0,1.0)+ where e err = error (show err)+ jac a b = abs $ (a!0) * (b!1) - (a!1) * (b!0)++{-|+ Integrates using a simple grid-sum.+-}+bzIntegral :: (Wavevector -> Parameterized Double) -> Parameterized Double+bzIntegral f = do kGrid <- meshBZ+ (Spacing dk) <- getMesh+ let df = liftM ((*) dk . (*) dk) . f+ liftM sum $! traverse df kGrid++{-|+ Returns a mesh of points within the /n/-parallelepiped subtended by the+ reciprocal lattice vectors. Covers the entire brillouin zone, though not+ in the shape you might expect, and not in a way that's pretty for graphing.+-}+meshBZ :: Parameterized (Grid Wavevector)+meshBZ = do bs <- recipPrimitiveLattice+ (Spacing dk) <- getMesh+ let frame = map (makeAxis dk) bs + return . foldl' (flip (uncurry M.insert)) M.empty + $ populate frame++{-|+ As 'meshBZ', but with parallelepipeds extending in each quadrant, octant...+ /n/-ant of the reciprocal lattice basis. Should cover more than the entire+ first Brillouin zone.+-}+meshBigBZ :: Parameterized (Grid Wavevector)+meshBigBZ = do bs <- recipPrimitiveLattice+ (Spacing dk) <- getMesh+ let frame = map (makeDoubleAxis dk) bs + return . foldl' (flip (uncurry M.insert)) M.empty + $ populate frame++--++makeAxis :: Double -> Wavevector -> [Wavevector]+makeAxis dk b = pos ++ [b]+ where pos = [ scale n bdk + | n <- takeWhile (\m -> (m*dk) < (norm_2 b)) [0..]]+ bdk = scale (dk / (norm_2 b)) b++makeDoubleAxis :: Double -> Wavevector -> [Wavevector]+makeDoubleAxis dk b = (reverse neg) ++ pos ++ [b]+ where pos = [ scale n bdk + | n <- takeWhile (\m -> (m*dk) < (norm_2 b)) [0..]]+ neg = [ scale n bdk + | n <- takeWhile (\m -> (abs $ m*dk) < (norm_2 b)) nins]+ bdk = scale (dk / (norm_2 b)) b+ nins = map negate [1..]++populate :: [[Wavevector]] -> [(GridIndex, Wavevector)]+populate [] = []+populate (ks:[]) = zipWith (\n k -> (GID [n],k)) [0..] ks+populate (ks:kss) = concatMap (\(GID i,k') -> [ (GID (j:i), k + k') + | (j,k) <- zip [0..] ks]) + $ populate kss
+ TBit/Hamiltonian/Builder/Decompactification.hs view
@@ -0,0 +1,74 @@+module TBit.Hamiltonian.Builder.Decompactification ( decompactify+ ) where++import TBit.Types+import TBit.Hamiltonian.Builder.PrimitiveLattice+import Data.Graph.Inductive+import Control.Monad.State (modify)++{-|+ Perform so-called "truncated decompactification" on a 'CellGraph'.++ Since the neighbor-data is stored in a unit-cell-level graph, it's in a+ sense "compact", i.e. it's a local periodic structure instead of an+ extended (to infinity) structure. Truncated decompactification sends the+ periodic structure (on T^2, roughly speaking) back to something of+ infinite extent (i.e. the integers), but then truncates the result to keep+ only a finite subset (i.e. the ribbon-width).++ It may not be clear /a priori/ how to choose the edge you want to+ 'decompactify' on to get the desired edge configuration; for honeycomb, + you can show on paper that decompactifying on a single graph edge (there+ are three, corresponding to the three nearest neighbors of a site) gives+ you zig-zag edge, while decompactifying on two graph edges gives you+ an armchair configuration. The square lattice is even more straightforward.+-}+decompactify :: Int + -> LEdge Displacement + -> CellGraph + -> Parameterizable CellGraph+decompactify n (v1,v2,d) gr = + do recordDecomEdges new+ setPrimLattice ret + return ret+ where g = replicateG n gr+ --les'= (replicateE n (noNodes gr)) (v1,v2,d)+ les'= concatMap (replicateE n (noNodes gr)) + $ filter (\(u,v,e) -> e == d) $ labEdges gr+ les = les' +++ map (\(x,y,r) -> (y,x,negate r)) les'+ new = filter (\(u,v,d) -> gelem u g && gelem v g) + $ map boost les+ ret = insEdges new $ delLEdges les g+ boost (u,v,e) = if e == d+ then (u, v + noNodes gr, e)+ else (u + noNodes gr, v, e)+{-+decompactify' :: Int + -> [LEdge Displacement]+ -> CellGraph + -> Parameterizable CellGraph+decompactify' n es gr = + do recordDecomEdges new+ setPrimLattice ret + return ret+ where g = replicateG n gr+ --les'= (replicateE n (noNodes gr)) (v1,v2,d)+ les'= concatMap (replicateE n (noNodes gr)) + $ filter (\(u,v,e) -> e == d) $ labEdges gr+ les = les' +++ map (\(x,y,r) -> (y,x,negate r)) les'+ new = filter (\(u,v,d) -> gelem u g && gelem v g) + $ map boost les+ ret = insEdges new $ delLEdges les g+ boost (u,v,e) = if e == d+ then (u, v + noNodes gr, e)+ else (u + noNodes gr, v, e)+-}+recordDecomEdges :: [LEdge Displacement] -> Parameterizable ()+recordDecomEdges es = modify (\ps -> ps { decomData = es ++ decomData ps })++(+++) :: [a] -> [a] -> [a]+(+++) [] [] = []+(+++) as [] = as+(+++) [] as = as+(+++) as bs = head as : (bs +++ tail as)+
+ TBit/Hamiltonian/Builder/Examples.hs view
@@ -0,0 +1,93 @@+module TBit.Hamiltonian.Builder.Examples where++import Data.Graph.Inductive.Graph+import Numeric.LinearAlgebra.HMatrix+import TBit.Types+import Data.Monoid++ring :: Int -> CellGraph+ring n = hermitize $ mkGraph vs es+ where vs = [ (j, ScalarSite j) | j <- [1..n] ]+ es = (n, 1, unit 1 - unit n) + : [ (j, j+1, unit (j+1) - unit j) | j <- [1..pred n]]+ m = fromIntegral n++ unit :: Int -> Vector Double+ unit j' = let j = fromIntegral j'+ in vector [ cos $ 2*pi*j/m , sin $ 2*pi*j/m ]++squareLattice :: CellGraph+squareLattice = mkGraph vs es+ where vs = [ (1, ScalarSite 1) ]+ es = [ (1, 1, vector [0,1])+ , (1, 1, vector [1,0])+ , (1, 1, vector [0,-1])+ , (1, 1, vector [-1,0]) ]++hexLattice :: CellGraph+hexLattice = hermitize $ mkGraph vs es+ where vs = [ (1, ScalarSite 1)+ , (2, ScalarSite 2)]+ es = [ (1, 2, vector [cos $ 2*pi/3, sin $ 2*pi/3])+ , (1, 2, vector [cos $ 4*pi/3, sin $ 4*pi/3]) + , (1, 2, vector [1,0]) ]++kagomeLattice :: CellGraph+kagomeLattice = hermitize $ mkGraph vs (map (fmap (scale 0.5)) es)+ where vs = [ (1, VectorSite 1 delta1)+ , (2, VectorSite 2 delta2)+ , (3, VectorSite 3 delta3)]+ es = [ (1, 2, vector [1,0])+ , (1, 2, negate $ vector [1,0])+ , (2, 3, vector [cos $ 4*pi/3, sin $ 4*pi/3]) + , (2, 3, negate $ vector [cos $ 4*pi/3, sin $ 4*pi/3]) + , (3, 1, vector [cos $ 2*pi/3, sin $ 2*pi/3])+ , (3, 1, negate $ vector [cos $ 2*pi/3, sin $ 2*pi/3]) ]+ delta1 = vector [ cos $ pi/6, sin . negate $ pi/6 ]+ delta2 = vector [ cos $ 7*pi/6 , sin $ 7*pi/6 ]+ delta3 = vector [ 0 , 1 ]+ ++instance Graph g => Monoid (g a b) where+ mempty = empty+ mappend g g' = mkGraph (labNodes g ++ labNodes g') + (labEdges g ++ labEdges g')++++kagomeRibbon :: Int -> CellGraph+kagomeRibbon = hermitize . kagomeRibbon'++kagomeRibbon' :: Int -> CellGraph+kagomeRibbon' 0 = hermitize $ mkGraph vs es+ where vs = [ (1, VectorSite 1 delta1)+ , (2, VectorSite 2 delta2) ]+ es = [ (1, 2, vector [0.5, 0]) ]+ delta1 = vector [ cos $ pi/6, sin . negate $ pi/6 ]+ delta2 = vector [ cos $ 7*pi/6 , sin $ 7*pi/6 ]+kagomeRibbon' n = (mkGraph vs es) `mappend` (boost $ kagomeRibbon (pred n))+ where vs = [ (1, VectorSite 1 delta1)+ , (2, VectorSite 2 delta2)+ , (3, VectorSite 3 delta3) ]+ es = [ (1, 2, e12)+ , (2, 3, e23)+ , (3, 1, e31)+ , (2, 1, e12)+ , (3, 4, negate e31)+ , (3, 5, e23)]+ e12 = vector [0.5, 0]+ e23 = vector [cos $ 4*pi/3, sin $ 4*pi/3]+ e31 = vector [cos $ 2*pi/3, sin $ 2*pi/3]+ delta1 = vector [ cos $ pi/6, sin . negate $ pi/6 ]+ delta2 = vector [ cos $ 7*pi/6 , sin $ 7*pi/6 ]+ delta3 = vector [ 0 , 1 ]+ boost g = let (e,v) = ( labEdges g , labNodes g )+ in mkGraph ( map (\(nd, VectorSite _ dt) + -> (nd + 3, VectorSite (nd + 3) dt)) v )+ ( map (\(src, tgt, dlt)+ -> (src + 3, tgt + 3, dlt)) e )++++hermitize :: CellGraph -> CellGraph+hermitize g = insEdges (map (\(x,y,r) -> (y,x,negate r)) $ labEdges g) g
+ TBit/Hamiltonian/Builder/Matrification.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE NoImplicitPrelude #-}+module TBit.Hamiltonian.Builder.Matrification (mixedSum, toMatrix) where++import Prelude.Listless+import Data.List.Stream+import Numeric.LinearAlgebra.HMatrix+import Data.Graph.Inductive+import Data.Maybe++mixedAdd :: Maybe (Matrix (Complex Double))+ -> Maybe (Matrix (Complex Double))+ -> Maybe (Matrix (Complex Double))+mixedAdd (Just a) (Just b) | dim a == dim b = Just $ a + b+ | dim a == 2 * dim b = Just $ a + expand b+ | dim b == 2 * dim a = Just $ expand a + b+ | otherwise = Nothing+ where dim = fst . size+ expand m = m `kronecker` ident 2++mixedAdd _ _ = Nothing++{-|+ Take a list of matrices, some of which may differ from the others in+ dimensionality by a factor of two, and maybe return the sum of these+ matrices with appropriate (right) kronecker products taken to make+ the summation well-formed. Useful for combining matrices written in+ forms with and without spin indices.+-}+mixedSum :: [Matrix (Complex Double)] -> Maybe (Matrix (Complex Double))+mixedSum = foldl1 mixedAdd . map Just++{-|+ Send a tight-binding graph model to the corresponding Hamiltonian+ matrix.+-}+toMatrix :: Gr (Matrix (Complex Double)) (Matrix (Complex Double)) + -> Matrix (Complex Double)+toMatrix gr = fromJust $ mixedSum [offd, ond]+ where es = map (\(i,j,k) -> ((i,j),k)) $ labEdges gr+ n = noNodes gr+ ns = labNodes gr+ d i j = if i==j then 1 else 0+ --offd = fromBlocks [[ fromMaybe 0 $ lookup (i,j) es | i <- [1..n]]+ -- | j <- [1..n]]+ offd = fromBlocks [[ sum $ map snd + $ filter ((==) (i,j) . fst) es | i <- [1..n]]+ | j <- [1..n]]+ ond = fromBlocks [[ (*) (d i j) + $ sum $ map snd $ filter ((==) j . fst) ns+ | i <- [1..n]]+ | j <- [1..n]]
+ TBit/Hamiltonian/Builder/PrimitiveLattice.hs view
@@ -0,0 +1,105 @@+module TBit.Hamiltonian.Builder.PrimitiveLattice ( primLattice+ , setPrimLattice + , delLEdges + , replicateE+ , replicateG+ ) where++import TBit.Types+import Numeric.LinearAlgebra.HMatrix+import Data.Graph.Inductive+import Data.Graph.Inductive.Query.SP (spLength)+import Data.Graph.Inductive.Query.DFS (components)+import Control.Monad.State+import Data.List.Stream (nub, nubBy, sortBy)++nthNeighborsTo :: CellGraph -> Int -> Node -> [(Node,Displacement)]+nthNeighborsTo g 0 j = [(j,vector [0,0])]+nthNeighborsTo g n j = concatMap (branch g) $ nthNeighborsTo g (pred n) j++branch :: CellGraph -> (Node, Displacement) -> [(Node, Displacement)]+branch g (i,r) = map (\(_,m,d) -> (m,r+d)) $ out g i+++{-|+ Sets the 'TBit.Types.latticeData' field of the 'TBit.Types.Parameters'+ state according to a 'primLattice'.+-}+setPrimLattice :: CellGraph -> Parameterizable CellGraph+setPrimLattice g = do tds <- gets decomData+ modify (\ps -> ps {latticeData = primLattice tds g})+ return g++{-|+ Determine a primitive lattice for a given 'TBit.Types.CellGraph'.+ + This is currently accomplished by determining the diameter of the thegraph,+ collecting all non-zero displacements from an arbitary site to itself+ which are within a diameter's worth of NN hopping, and returning a+ maximal linearly independent subset of these with a preference for+ vectors with smaller L2 norms.++ For finite nanoribbons such as those generated by+ 'TBit.Hamiltonian.Builder.Decompactification.decompactify', the graph+ diameter can be quite large. This means that finding the primitive+ lattice vectors as described above can become unreasonably slow. In the+ future, we will label decompactified lattice vectors in the CellGraph+ so that 'graphDiameter' will understand not to count them.+-}+primLattice :: [LEdge Displacement] -> CellGraph -> Lattice+primLattice tds gr = take (rank' $ fromColumns $ bulkCycles tds gr site) + $ bulkCycles tds gr site+ where site = fst . head . labNodes $ gr+ rank' m = if fst (size m) == 0+ then 0+ else rank m++bulkCycles :: [LEdge Displacement] -> CellGraph -> Node -> Lattice+bulkCycles tds gr site = nub+ $ nubBy (\u v -> rank (fromColumns [u,v]) /= 2)+ $ sortBy (\a b -> compare (norm_2 a) (norm_2 b)) + $ filter ((<) 0.0001 . norm_2)+ $ map snd twins+ where hops = nthNeighborsTo gr (succ $ graphDiameter g) site + ges = filter (\(u,v,_) -> (u `elem` gvs)+ && (v `elem` gvs)) $ labEdges gr+ gvs = head $ components $ delLEdges tds gr+ g :: Gr Int Int+ g = emap (const 1) $ mkGraph (map (\x -> (x,x)) gvs) ges+ twins = filter next hops+ next (k,d) = (site == k) && (any (/= 0) $ toList d)++graphDiameter :: Graph gr => gr a Int -> Int+graphDiameter g = maximum $ [ spLength i j g | i <- [1..noNodes g]+ , j <- [1..noNodes g]+ , i <= j ]+++-- | Delete a list of 'Data.Graph.Inductive.LEdge's from a graph. +delLEdges :: (Eq b, DynGraph gr) => [LEdge b] -> gr a b -> gr a b+delLEdges es g = foldr delLEdge g es++-- | Replicate a 'Data.Graph.Inductive.LEdge' n times, each time+-- increasing the in and out nodes by m. (n is the first argument,+-- m the second, somewhat stupidly.)+replicateE :: Int -> Int -> LEdge Displacement -> [LEdge Displacement]+replicateE n m (v1,v2,e) = [ (v1 + j*m, v2 + j*m, e) | j <- [0..pred n]]++-- | Replicate a 'TBit.Types.CellGraph' n times, each time+-- increasing the in and out nodes by the number of nodes+-- in the 'TBit.Types.CellGraph'.+replicateG :: Int -> CellGraph -> CellGraph+replicateG n gr = mkGraph vss ess+ where m = noNodes gr+ vss = concatMap (\(j,vs) + -> map (\(v,a) + -> (v + j*m, a { num = num a + j*m })) vs)+ $ zip [0..] + $ map labNodes gs+ ess = concatMap (\(j,es) + -> map (\(v1,v2,r) + -> (v1 + j*m, v2 + j*m, r)) es)+ $ zip [0..] + $ map labEdges gs+ gs = replicate n gr+
+ TBit/Hamiltonian/Builder/Terms.hs view
@@ -0,0 +1,138 @@+{-# LANGUAGE NoImplicitPrelude #-}+module TBit.Hamiltonian.Builder.Terms ( neighborTerm+ , onsiteTerm+ , parityStaggeredTerm+ , localMoments+ , kagomeSOC+ , rashbaZ+ ) where++import Prelude.Listless+import Data.List.Stream+import TBit.Types+import TBit.Parameterization+import Control.Monad.State+import Data.Graph.Inductive+import Data.Maybe+import Numeric.LinearAlgebra.HMatrix++-- | Add nearest-neighbor hopping to a lattice model.+neighborTerm :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+neighborTerm p g = do t <- getScalar p+ r <- restrictor+ return (\k -> nmap (scalar . const 0) + $ emap (scalar . \d -> t * cis ((r d) `dot` (r k))) g)++-- | Add an onsite energy term to a lattice model.+onsiteTerm :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+onsiteTerm p g = do t <- getScalar p+ return (\_ -> nmap (scalar . const t) + $ emap (scalar . const 0) g)++-- | Add an staggered onsite term to a lattice model. Works based on the+-- integer parity of graph nodes, making it model-detail-dependent.+parityStaggeredTerm :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+parityStaggeredTerm p g = do t <- getScalar p+ return (\_ -> nmap (scalar . \e -> t * (-1.0)^(num e + 1)) + $ emap (scalar . const 0) g)+{-|+ Produces a representation of local magnetic moments given by+ site-wise 'TBit.Types.VectorSite' data. Fails clumsily if+ applied to Scalar sites.+-}+localMoments :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+localMoments p g = do t <- getScalar p+ return (\_ -> nmap (scale t . exch . mom)+ $ emap (const zeroMtx) g)+ where zeroMtx = scale 0 (ident 2)+ exch = sum + . zipWith (flip scale) [sigmaX, sigmaY, sigmaZ]+ . map (:+ 0.0)+ . toList++{-| Produces a Rashba spin-orbit coupling term for an E-field applied along+ the /z/ direction. -}+rashbaZ :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+rashbaZ p g = do t <- getScalar p+ r <- restrictor+ return $ \k -> mtx t k r+ where nng = nmap (\u -> map ((,) (num u))+ . filter (\t -> snd t /= vector [0,0]) + $ nthNeighborsTo g 1 (num u)) g+ rsh t k r = nmap + ( map (\(u,(v,d)) -> (,,) u v+ . scale (iC * t)+ . scale (cis $ r k `dot` r d)+ $ scale (c $ d!0) sigmaY+ - scale (c $ d!1) sigmaX)) nng+ mtx t k r = let ns = labNodes $ rsh t k r+ es = concatMap snd ns+ in nmap (const 0) $ mkGraph ns es+ c z = z :+ 0.0+ ++{-|+ Produces a spin-orbit interaction in the style of that given by+ Hua et al in PRL /112/, 017205 (2014). It should probably only+ be used with a Kagomé lattice 'TBit.Types.CellGraph'. It works by+ looking at nearest neighbor pairs (/i/,/j/) and then looking up the+ 'TBit.Types.VectorSite' moment of site /k/; /k/ is computed as+ + > k = let i' = succ $ (i - 1) `mod` 3+ > j' = succ $ (j - 1) `mod` 3+ > in head $ [1,2,3] \\ [i',j']++ (Recall that /i/ and /j/ are indexed from 1 as nodes.) Clearly, + this function is not safe unless it's applied to the correct lattice.++ Once 'TBit.Types.mom' is computed for 'TBit.Types.VectorSite' /k/, it+ is coupled to the Pauli matrix tensor as expected. The parity ν+ is chosen by asking whether succ /i/ == /j/ mod 3.+-}+kagomeSOC :: String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)+kagomeSOC p g = do t <- getScalar p+ r <- restrictor+ return $ (\k -> mtx t k r)+ -- The nng graph is just like g, but nodes are replaced with+ -- lists of nearest neighbors and their displacements+ where nng = nmap (\u -> map ((,) (num u))+ . filter (\t -> snd t /= 0) + $ nthNeighborsTo g 1 (num u)) g+ moms t k r = nmap + ( map (\(u,(v,d)) -> (,,) u v+ . scale (iC * t)+ . scale (cis $ (r k) `dot` (r d))+ $! exch (nab u v))) nng+ mtx t k r = let ns = labNodes $! moms t k r+ es = concatMap snd ns+ in nmap (const 0) $! mkGraph ns es+ nab i j = let i' = succ $ pred i `mod` 3+ j' = succ $ pred j `mod` 3+ in neg i j . mom . fromJust . lab g . head $! [1,2,3] \\ [i',j']+ neg i j = let i' = succ $ pred i `mod` 3+ j' = succ $ pred j `mod` 3+ in if (succ i' `mod` 3) == (j' `mod` 3)+ then scale 1.0+ else scale (-1.0)+ exch = sum + . zipWith (flip scale) [sigmaX, sigmaY, sigmaZ]+ . map (:+ 0.0)+ . toList+ + ++nthNeighborsTo :: CellGraph -> Int -> Node -> [(Node,Displacement)]+nthNeighborsTo _ 0 j = [(j,vector [0,0])]+nthNeighborsTo g n j = concatMap (branch g) $! nthNeighborsTo g (pred n) j++branch :: CellGraph -> (Node, Displacement) -> [(Node, Displacement)]+branch g (i,r) = map (\(_,m,d) -> (m,r+d)) $ out g i++restrictor :: Parameterized (Vector Double -> Vector Double)+restrictor = do (l:ls) <- gets latticeData+ case compare (size l) (rank . fromColumns $ l:ls)+ of EQ -> return id+ GT -> return $ (#>) (fromRows (map nrml $ l:ls))+ LT -> error "primitive lattice overdetermines space"+ where nrml v = scale (1.0/(norm_2 v)) v+
+ TBit/Hamiltonian/Eigenstates.hs view
@@ -0,0 +1,41 @@+module TBit.Hamiltonian.Eigenstates ( eigenstates+ , eigenbras+ , eigenkets+ , eigenenergies+ , eigensystem ) where++import TBit.Types+import Data.List.Stream (sortBy, zip, sort, map)+import Prelude hiding (zip, map)+import Numeric.LinearAlgebra.HMatrix++-- | Returns a list of eigenvectors sorted by eigenvalue. The lowest+-- energy state is the first element of the returned list.+eigenstates :: Hamiltonian -> Wavevector -> Parameterized [Eigenstate]+eigenstates h k = return . map snd $ eigensystem' (h k)++-- |As 'eigenstates', but converts each vector to a column matrix+-- for convenience in certain caclulations.+eigenkets :: Hamiltonian -> Wavevector -> Parameterized [Eigenket]+eigenkets h k = eigenstates h k >>= (return . map asColumn)++-- |As 'eigenkets', takes the conjugate transpose of each ket.+eigenbras :: Hamiltonian -> Wavevector -> Parameterized [Eigenbra]+eigenbras h k = eigenkets h k >>= (return . map tr)++-- |Returns a list of eigenvalues, sorted in ascending order.+eigenenergies :: Hamiltonian -> Wavevector -> Parameterized [Energy]+eigenenergies h = return . energies . h++energies :: Matrix ℂ -> [Energy]+energies = sort . toList . eigenvaluesSH'++eigensystem' :: Matrix ℂ -> [(Energy,Eigenstate)]+eigensystem' = sortBy (\t1 t2 -> compare (fst t1) (fst t2))+ . (\(vals,vecs) -> zip (toList vals) (toColumns vecs))+ . eigSH'++-- |Returns the full eigensystem, sorted by energy. Equivalent to zipping+-- the results of 'eigenenergies' and 'eigenstates'.+eigensystem :: Hamiltonian -> Wavevector -> Parameterized [(Energy,Eigenstate)]+eigensystem h k = return $ eigensystem' (h k)
+ TBit/Magnetic/OrbitalMagnetization.hs view
@@ -0,0 +1,144 @@+module TBit.Magnetic.OrbitalMagnetization ( orbMag+ , dichroism+ , dichroicIntegrand+ , intrinsicOM+ , integrandMk+ , bandIntrinsicOM+ ) where++import TBit.Types+import TBit.Framework+import TBit.Hamiltonian.Eigenstates+import TBit.Numerical.Derivative+import TBit.Parameterization++import Data.List.Stream (delete)+import Control.Monad.Except+import Prelude hiding (mapM, sequence)+import Numeric.LinearAlgebra.HMatrix++-- |Returns the total orbital magnetization due to filling the first+-- n bands.+orbMag :: Filling -> Hamiltonian -> Parameterized Magnetization+orbMag n h = bzIntegral (integrandOM n h)++-- |Returns the gauge-invariant self-rotational orbital magnetization, i.e.+-- the circular dichroism, of the occupied bands, accomplished by integrating+-- 'dichroicIntegrand'.+dichroism :: Filling -> Hamiltonian -> Parameterized Magnetization+dichroism n h = bzIntegral (dichroicIntegrand n h)++-- |Returns the intrinsic orbital magnetization of the n/th/ band, namely+-- the integral of _m_(k) from (Xiao et al., 2005).+bandIntrinsicOM :: BandIndex -> Hamiltonian -> Parameterized Magnetization+bandIntrinsicOM n h = bzIntegral (integrandMk n h)++-- |Sums 'bandIntrinsicOM' over the first n bands.+intrinsicOM :: Filling -> Hamiltonian -> Parameterized Magnetization+intrinsicOM n h = liftM sum $ sequence [ bandIntrinsicOM j h | j <- [0 .. pred n] ]++-- |Essentially equation (12) from PRB _77_, 054438 (2008) with alpha, beta set+-- to x, y respectively so as to give the z-component of the (expectation+-- value of) the circular dichroism pseudovector. This form takes a double sum over+-- occupied and then unoccupied states, which improves over the implementation of+-- 'integrandMk' by allowing for intra-(un)occupied band degeneracies as long as +-- there is still an electronic band gap.+-- +-- As is consistent with our API for these types of functions, the +-- 'TBit.Types.Filling' argument should be a positive integer counting the+-- number of filled bands.+dichroicIntegrand :: Filling -> Hamiltonian -> Wavevector -> Parameterized Magnetization+dichroicIntegrand b h k = do ket <- eigenkets h k+ bra <- eigenbras h k+ eng <- eigenenergies h k+ (Spacing s) <- getMesh+ return . negate . imagPart . sum+ $ zipWith (/)+ [ num $ (bra!!n) <> hx s <> (ket!!m)+ <> (bra!!m) <> hy s <> (ket!!n)+ | m <- occ , n <- unocc ]+ [ (eng!!m - eng!!n) :+ 0.0+ | m <- occ , n <- unocc ]+ where num m = m ! 0 ! 0+ occ = [0 .. pred b]+ unocc = [b .. pred dim]+ dim = fst $ size $ h k+ twice = (*) 2.0+ kx = k ! 0+ ky = k ! 1+ hx s = diff (s :+ 0.0) (\x -> h $ vector [realPart x, ky]) (kx :+ 0.0)+ hy s = diff (s :+ 0.0) (\y -> h $ vector [kx, realPart y]) (ky :+ 0.0)++ ++-- |Gives the gauge-invariant self-rotational orbital magnetism, which is+-- proportional to the circular dichroism for a particular band index.+-- This is m(k) in Xiao et al's semiclassical approach (hence the name).+integrandMk :: BandIndex -> Hamiltonian -> Wavevector -> Parameterized Magnetization+integrandMk n h k = do ket <- eigenkets h k+ bra <- eigenbras h k+ eng <- eigenenergies h k+ (Spacing s) <- getMesh+ if (close eng n (pred n)) || (close eng n (succ n))+ then throwError $ UndefinedError (degenErr eng)+ else return + $ imagPart . sum+ $ zipWith (/)+ [ num $! (bra!!m) <> hx s <> (ket!!n)+ <> (bra!!n) <> hy s <> (ket!!m)+ | m <- delete n [0 .. pred dim]]+ [ (eng!!n - eng!!m) :+ 0.0+ | m <- delete n [0 .. pred dim]]+ where num m = m ! 0 ! 0+ dim = fst $ size $ h k+ kx = k ! 0+ ky = k ! 1+ eps = 0.000000001+ hx s = diff (s :+ 0.0) (\x -> h $ vector [realPart x, ky]) (kx :+ 0.0)+ hy s = diff (s :+ 0.0) (\y -> h $ vector [kx, realPart y]) (ky :+ 0.0)+ close _ (-1) _ = False+ close _ _ (-1) = False+ close eng i j = (abs $ eng!!i - eng!!j) < eps+ degenErr eng = case (close eng n (pred n), close eng n (succ n))+ of (True, _) -> errMsg $ pred n+ (_, True) -> errMsg $ succ n+ otherwise -> "This error was in error."+ errMsg j = "Cannot integrate Berry curvature due to a gap "+ ++"closing between bands "++(show $ succ j)++" and "+ ++(show $ succ n)++"."++-- |The full orbital magnetization for a particular wavevector; integrating over+-- the Brillouin zone would give the actual OM due to a particular band filling.+integrandOM :: Filling -> Hamiltonian -> Wavevector -> Parameterized Magnetization+integrandOM n h k = do ket <- eigenkets h k+ bra <- eigenbras h k+ eng <- eigenenergies h k+ (Spacing s) <- getMesh+ if any (< eps) (diffs eng)+ then throwError $ UndefinedError (degenErr eng)+ else return + $ imagPart . sum+ $ zipWith (/)+ [ (*) ((eng!!m + eng!!l) :+ 0.0) .+ num $! (bra!!m) <> hx s <> (ket!!l)+ <> (bra!!l) <> hy s <> (ket!!m)+ | m <- occ , l <- unocc ]+ [ (eng!!m - eng!!l)^2 :+ 0.0+ | m <- occ , l <- unocc ]+ where num m = m ! 0 ! 0+ occ = [0 .. pred n]+ unocc = [n .. pred dim]+ dim = fst $ size $ h k+ diffs eng = [ abs $ (eng!!m) - (eng!!l)+ | m <- occ, l <- unocc ]+ kx = k ! 0+ ky = k ! 1+ eps = 0.000000001+ hx s = diff (s :+ 0.0) (\x -> h $ vector [realPart x, ky]) (kx :+ 0.0)+ hy s = diff (s :+ 0.0) (\y -> h $ vector [kx, realPart y]) (ky :+ 0.0)+ degenErr eng = let (j,l,_) = head $ filter (\(_,_,z) -> (z < eps))+ $ [ (m,p,abs $ eng!!m - eng!!p)+ | m <- occ, p <- unocc ]+ in "Cannot integrate Berry curvature due to a gap "+ ++"closing between bands "++(show $ succ j)++" and "+ ++(show $ succ l)++"."
+ TBit/Numerical/Derivative.hs view
@@ -0,0 +1,25 @@+module TBit.Numerical.Derivative (diff, diff4) where++import Numeric.LinearAlgebra.HMatrix++-- |Takes the element-wise first derivative of the given vector/matrix-valued+-- function of a single variable. The first argument gives the spacing epsilon+-- as used in the difference formula, /i.e./ that thing which limit is taken+-- to zero. The second argument is the matrix function, and the third+-- argument is the point at which to evaluate the derivative. +--+-- The derivative is taken using the finite difference method to second order.+diff :: (Num (c e), Container c e) => e -> (e -> c e) -> e -> c e+diff e m x = sum+ [ scale (2.0/3.0/e) (m (x+e) - m (x-e))+ , scale (1.0/12.0/e) (m (x-2*e) - m (x+2*e)) ]++-- |As 'diff', but uses fourth order finite difference method. This means that+-- the matrix argument will need to be evaluated at twice as many points, which+-- increases the runtime twofold until we can implement a full-gridded calculation.+diff4 :: (Num (c e), Container c e) => e -> (e -> c e) -> e -> c e+diff4 e m x = sum+ [ scale (4.0/5.0/e) (m (x+e) - m (x-e))+ , scale (1.0/5.0/e) (m (x-2*e) - m (x+2*e))+ , scale (4.0/105.0/e) (m (x+3*e) - m (x-3*e))+ , scale (1.0/280.0/e) (m (x-4*e) - m (x+4*e))]
+ TBit/Numerical/Integration.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}++module TBit.Numerical.Integration (integrate) where++import Numeric.Integration.TanhSinh++type Res = (Double, Double)+type Bound = (Double, Double)++ttrap :: (Double -> Double) -> Double -> Double -> Res+ttrap f xmin xmax = (ans, err)+ where+ res = absolute 1e-3 $ parTrap f xmin xmax+ ans = result res+ err = errorEstimate res++class Integrable r where+ type Bounds r :: *+ integrate' :: r -> Bounds r -> Res++instance Integrable Double where+ type Bounds Double = ()+ integrate' x _ = (x, 0)++instance Integrable r => Integrable (Double -> r) where+ type Bounds (Double -> r) = (Bound, Bounds r)+ integrate' f ((xmin, xmax), args) = ttrap g xmin xmax+ where+ g x = fst $ integrate' (f x) args++class CurryBounds bs where+ type Curried bs a :: *+ curryBounds :: (bs -> a) -> Curried bs a++instance CurryBounds () where+ type Curried () a = a+ curryBounds f = f ()++instance CurryBounds bs => CurryBounds (b, bs) where+ type Curried (b, bs) a = b -> Curried bs a+ curryBounds f = \x -> curryBounds (\xs -> f (x, xs))++-- |Integrate a function of n variables by giving the corresponding n integration+-- domains, /i.e./+--+-- > let f x y z = x^2 + log (y-z)+-- > integrate f (0,1) (3,4) (1,2)+--+-- This code was borrowed from the excellent answer at http://stackoverflow.com/questions/23703360/using-numeric-integration-tanhsinh-for-n-dimensional-integration.+-- +-- The integration uses the Tanh-Sinh quadrature method and relies on Kmett's integration libary.+integrate :: (Integrable r, CurryBounds (Bounds r)) => r -> Curried (Bounds r) Res+integrate = curryBounds . integrate'
+ TBit/Parameterization.hs view
@@ -0,0 +1,88 @@+module TBit.Parameterization ( loadParams+ , getScalar+ , getVector+ , getMesh+ , crunch+ , primitiveLattice+ , recipPrimitiveLattice) where++import TBit.Types+import Prelude hiding (map, head)++import Control.Monad.State+import Control.Monad.Except ++import qualified Data.Map as M+import Data.Complex+import Data.List.Stream (map, head)+import Data.Maybe++import Numeric.LinearAlgebra.HMatrix (toColumns, fromRows, col, linearSolve, Vector)++-- |Given a list of ('Prelude.String', a) pairs, return a mapping from the+-- string to the value. Such a map is suitable for setting the 'TBit.Types.scalarParams'+-- and 'TBit.Types.vectorParams' records of the 'TBit.Types.Parameters' type.+loadParams :: [(String, a)] -> M.Map String a+loadParams = foldl (flip (uncurry M.insert)) M.empty++-- |Given a 'Prelude.String', retrieve that value from the parameterization monad.+-- If no such scalar has been uploaded to the 'TBit.Types.Parameters' type, 'getScalar'+-- with throw an error in 'TBit.Types.Parameterized'\'s 'Control.Monad.Except.ExceptT' monad.+getScalar :: String -> Parameterized (Complex Double)+getScalar s = do p <- gets (M.lookup s . scalarParams)+ if isJust p+ then return $ fromJust p+ else throwError + $ UnknownParameter s++realScalar :: String -> Parameterized Double+realScalar = liftM realPart . getScalar++complexScalar :: String -> Parameterized (Complex Double)+complexScalar = getScalar++-- |Given a 'Prelude.String', retrieve that value from the parameterization monad.+-- If no such scalar has been uploaded to the 'TBit.Types.Parameters' type, 'getScalar'+-- with throw an error in 'TBit.Types.Parameterized'\'s 'Control.Monad.Except.ExceptT' monad.+getVector :: String -> Parameterized (Vector (Complex Double))+getVector s = do p <- gets (M.lookup s . vectorParams)+ if isJust p+ then return $ fromJust p+ else throwError + $ UnknownParameter s++-- |Return the spacing information for purposes of gridding the Brillouin zone.+getMesh :: Parameterized Meshing+getMesh = gets meshingData++-- |Compute a 'TBit.Types.Parameterized' quantity given a set of input parameters.+crunch :: Parameterized a -> Parameters -> Either TBError a+crunch pmz ps = fst $ runState (runExceptT pmz) ps++-- |Return a list of primitive lattice vectors.+primitiveLattice :: Parameterized Lattice+primitiveLattice = gets latticeData++-- |Return a list of reciprocal primitive lattice vectors. They correspond in order+-- to the return values of 'primitiveLattice', in the sense that:+-- +-- > do as <- primitiveLattice+-- > bs <- primitiveLattice+-- > return $ zipWith `dot` as bs+--+-- will return the list+--+-- > replicate dim (2*pi)+--+-- up to numerical precision.+recipPrimitiveLattice :: Parameterized Lattice+recipPrimitiveLattice = do lat <- primitiveLattice+ let l = length lat+ let bs = [ linearSolve (fromRows lat) + (col $ mtx l n)+ | n <- [1..l]]+ if any isNothing bs+ then throwError SingularLatticeError+ else (return . map (head . toColumns) . catMaybes) bs+ where mtx l n = [ 2.0 * pi * kd n j | j <- [1..l] ]+ kd a b = if a == b then 1.0 else 0.0
+ TBit/Plots.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE NoImplicitPrelude #-}+module TBit.Plots (bzPlot, paramPlot, bandPlot) where++import TBit.Types+import TBit.Framework+import TBit.Parameterization+import TBit.Hamiltonian.Eigenstates (eigenenergies)++import Prelude.Listless+import Data.List.Stream+import Data.Map (adjust, elems)+import Control.Monad.State+import Control.Parallel.Strategies+import Numeric.LinearAlgebra.HMatrix+import qualified Data.Traversable as T++{-|+ Given a hamiltonian and a set of parameters, plot a real-valued function+ over the Brillouin zone. For a Berry curvature plot, /e.g./,+ + > bzPlot (kagomeAF, defaultParams) (bandCurvature 0)++ The output string is suitable for gnuplot, unless that string is+ reporting an error from the calculation.+-}+bzPlot :: (Parameterized Hamiltonian, Parameters) + -> (Hamiltonian -> Wavevector -> Parameterized Double) + -> String+bzPlot (h, ps) f = case crunch (h >>= computeBZPlot f) ps+ of Left err -> show err+ Right dat -> unlines+ . map (\(k, d) -> (unwords . map show . toList $ k)+ ++ " " ++ show d)+ $ dat++computeBZPlot :: (Hamiltonian -> Wavevector -> Parameterized Double)+ -> Hamiltonian + -> Parameterized [(Wavevector, Double)]+computeBZPlot f h = meshBZ + >>= T.mapM (\k -> T.sequence (k, f h k))+ >>= (return . elems)++{-|+ Given a hamiltonian and a set of parameters, plot a real-valued function+ over some range of a parameter. For a Chern number plot over the hopping+ parameter, /e.g./,+ + > paramPlot (kagomeAF, defaultParams) (chern 1) ("t", [1.0, 1.1 .. 5.0])++ The output string is suitable for gnuplot, unless that string is+ reporting an error from the calculation.+-}+paramPlot :: (Parameterized Hamiltonian, Parameters)+ -> (Hamiltonian -> Parameterized Double)+ -> (String, [Double])+ -> String+paramPlot (h, ps) f (s, ss) = case crunch (computePMPlot f s ss h) ps+ of Left err -> show err+ Right dat -> unlines+ . map (\(p, d) -> show p + ++ " " + ++ show d)+ $ dat+ ++computePMPlot :: (Hamiltonian -> Parameterized Double)+ -> String+ -> [Double]+ -> Parameterized Hamiltonian+ -> Parameterized [(Double, Double)]+computePMPlot f s ss h = liftM (zip ss) + $ mapM (\v -> modify (setScalar s v) >> (h >>= f)) $ ss+ where setScalar str val ps = ps { scalarParams = adjust (\_ -> val :+ 0.0) str + $ scalarParams ps }+++multiPlot :: (Parameterized Hamiltonian, Parameters)+ -> [(Hamiltonian -> Parameterized Double)]+ -> (String, [Double])+ -> String+multiPlot sys fs params = unlines $ map (\f -> paramPlot sys f params) fs++{-|+ Given a hamiltonian and a set of parameters, plot the bands over+ a path anchored by the provided wavevectors. For instance,+ + > gamma = vector [0.0 , 0.0]+ > point1 = vector [0.5 , 1.0]+ > point2 = vector [0.5 , -1.0]+ > bandPlot (kagomeAF, defaultParams) [gamma, point1, point2, gamma]++ The output string is suitable for gnuplot, unless that string is+ reporting an error from the calculation.+-}+bandPlot :: (Parameterized Hamiltonian, Parameters)+ -> [Wavevector]+ -> String+bandPlot (h,ps) ks = either show id $ crunch (bandData h ks) ps++bandData :: Parameterized Hamiltonian+ -> [Wavevector]+ -> Parameterized String+bandData h ks = do mtx <- h+ kPath ks >>= mapM (eigenenergies mtx)+ >>= (return . unlines+ . map (\(x, e) -> show x ++ " " ++ show e)+ . concatMap T.sequence + . zip [0..])
+ TBit/Systems/HoneycombLattice.hs view
@@ -0,0 +1,64 @@+module TBit.Systems.HoneycombLattice (defaultParams, kaneMele, parameters) where++import TBit.Types+import TBit.Parameterization +import Numeric.LinearAlgebra.HMatrix+import qualified Data.Map as M++-- |The default set of scalar parameters is:+-- +-- * "t" = 1, the hopping parameter+-- * "soc" = 1, the intrinsic spin-orbit coupling+-- * "r" = 1, the Rashba spin-orbit coupling+-- * "v" = 1, the staggered on-site energy+-- * "hz" = 1, the out-of-plane AF parameter+--+-- These can be changed by using 'parameters' to set+-- all of them explicitly.+defaultParams :: Parameters+defaultParams = parameters [ ("t" , 1.0 :+ 0.0)+ , ("soc", 1.0 :+ 0.0)+ , ("r" , 1.0 :+ 0.0)+ , ("v" , 1.0 :+ 0.0)+ , ("hz" , 1.0 :+ 0.0) ]++-- | Set the named parameters to their given complex values.+-- This function is used to implement 'defaultParams' as+-- +-- > defaultParams = parameters [ ("t" , 1.0 :+ 0.0)+-- > , ("soc", 1.0 :+ 0.0)+-- > , ("r" , 1.0 :+ 0.0)+-- > , ("v" , 1.0 :+ 0.0)+-- > , ("hz" , 1.0 :+ 0.0) ]+--+-- but you can use it to generate your own parameter list. For+-- more advanced manipulation, like setting the mesh size or the+-- primitive lattice vectors, you'll have to constuct a 'Types.Parameters'+-- type explicitly.+parameters :: [(String, Complex Double)] -> Parameters+parameters ps = Parameters { latticeData = [ vector [ 0.5, 0.5 * sqrt 3] + , vector [-0.5, 0.5 * sqrt 3] ] + , scalarParams = loadParams ps+ , vectorParams = M.empty + , meshingData = Spacing 0.1 }++-- | Exported directly from a Mathematica and hard-coded into this module.+-- Uses the exact conventions of Kane Mele, excep that the Gamma_15 term+-- also includes a constant value for AF order. This AF order is taken+-- in the large Hubbard U limit by construction.+kaneMele :: Parameterized Hamiltonian+kaneMele = do+ t <- getScalar "t"+ v <- getScalar "v"+ hz <- getScalar "hz"+ r <- getScalar "r"+ soc <- getScalar "soc"+ return $ \k -> let kx = (k!0) :+ 0.0+ ky = (k!1) :+ 0.0+ y = realPart ky+ in (4 >< 4) + [ hz + v - 4*soc*cos((sqrt(3)*ky)/2)*sin(kx/2) + 2*soc*sin(kx) , 0, t + (2*t*cos(kx/2))/cis(0.5*sqrt(3)*y), r*(negate iC + (iC*(cos(kx/2) + sqrt(3)*sin(kx/2)))/cis(0.5*sqrt(3)*y))+ , 0, -hz + v + 4*soc*cos((sqrt(3)*ky)/2)*sin(kx/2) - 2*soc*sin(kx), r*(negate iC + (iC*(cos(kx/2) - sqrt(3)*sin(kx/2)))/cis(0.5*sqrt(3)*y)), t + (2*t*cos(kx/2))/cis(0.5*sqrt(3)*y) + , t + 2*cis(0.5*sqrt(3)*y)*t*cos(kx/2), r*(iC - (negate iC)*cis(0.5*sqrt(3)*y)*(cos(kx/2) - sqrt(3)*sin(kx/2))), -hz - v + 4*soc*cos((sqrt(3)*ky)/2)*sin(kx/2) - 2*soc*sin(kx), 0+ , r*(iC - iC*cis(0.5*sqrt(3)*y)*(cos(kx/2) + sqrt(3)*sin(kx/2))), t + 2*cis(0.5*sqrt(3)*y)*t*cos(kx/2), 0, hz - v - 4*soc*cos((sqrt(3)*ky)/2)*sin(kx/2) + 2*soc*sin(kx) ]+
+ TBit/Systems/KagomeLattice.hs view
@@ -0,0 +1,132 @@+module TBit.Systems.KagomeLattice (defaultParams, parameters, kagomeAF) where++import Control.Monad (liftM2)+import Numeric.LinearAlgebra.HMatrix+import TBit.Types+import TBit.Parameterization+import Data.Foldable hiding (sum, toList)+import Data.Map (empty) ++-- |The default set of scalar parameters is:+-- +-- * "t" = 1, the hopping parameter+-- * "tSO" = 1, the intrinsic spin-orbit coupling+-- * "J" = 1, the Heisenberg exchange parameter+--+-- The default vector parameters are the d-orbital local moments+-- on the three sites. Each of them takes the form:+-- (-cos theta, -sin theta, 0) where theta is:+--+-- * "d0" : theta = pi/2+-- * "d1" : theta = pi/2 + 2pi/3+-- * "d2" : theta = pi/2 + 4pi/3+--+-- These can be changed by using 'parameters' to set+-- all of them explicitly.+defaultParams = parameters [ ("t" , 1.0 :+ 0.0)+ , ("tSO", 0.2 :+ 0.0)+ , ("J" , 1.7 :+ 0.0) ]+ [ ("d0" , n21 )+ , ("d1" , n02 )+ , ("d2" , n10 ) ]+ where n10 = negate $ fromList [ cos ang01 , sin ang01 , 0.0 ]+ n21 = negate $ fromList [ cos ang12 , sin ang12 , 0.0 ]+ n02 = negate $ fromList [ cos ang20 , sin ang20 , 0.0 ]+ ang01 = pi/2.0 + 4.0*pi/3.0+ ang12 = pi/2.0+ ang20 = pi/2.0 + 2.0*pi/3.0++-- | Set the named parameters to their given complex values.+-- This function is used to implement 'defaultParams' as+--+-- > defaultParams = parameters [ ("t" , 1.0 :+ 0.0)+-- > , ("tSO", 0.2 :+ 0.0)+-- > , ("J" , 1.7 :+ 0.0) ]+-- > [ ("d0" , n21 )+-- > , ("d1" , n02 )+-- > , ("d2" , n10 ) ]+-- > where n10 = negate $ fromList [ cos ang01 , sin ang01 , 0.0 ]+-- > n21 = negate $ fromList [ cos ang12 , sin ang12 , 0.0 ]+-- > n02 = negate $ fromList [ cos ang20 , sin ang20 , 0.0 ]+-- > ang01 = pi/2.0 + 4.0*pi/3.0+-- > ang12 = pi/2.0+-- > ang20 = pi/2.0 + 2.0*pi/3.0+--+-- but you can use it to generate your own parameter list. For+-- more advanced manipulation, like setting the mesh size or the+-- primitive lattice vectors, you'll have to constuct a 'Types.Parameters'+-- type explicitly.+parameters :: [(String, Complex Double)] -> [(String, Vector (Complex Double))] -> Parameters+parameters ps vs = Parameters { latticeData = [ vector [ 1.0, 0.0 ]+ , vector [cos (pi/3.0),sin (pi/3.0)]]+ , scalarParams = loadParams ps+ , vectorParams = loadParams vs+ , meshingData = Spacing 0.1 }+++-- |The kagomé hamiltonian provided here includes nearest-neighbor hopping,+-- noncollinear AF order due to localized d-orbital moments, and spin-orbit+-- coupling which breaks mirror symmetry.+kagomeAF :: Parameterized Hamiltonian+kagomeAF = do hop <- hopping+ af <- afOrder+ soc <- spinOrbit+ return $ \k -> hop k + af k + soc k++hopping :: Parameterized Hamiltonian+hopping = do t <- getScalar "t"+ lat <- primitiveLattice+ return $ \k -> flip kronecker (ident 2)+ . plusCT+ . scale t+ . (3 >< 3)+ $ [ 0 , phase01 k lat , 0+ , 0 , 0 , phase12 k lat+ , phase20 k lat , 0 , 0 ]+ where plusCT m = m + tr m+ phase01 k (_:a2:_) = (2 * cos(0.5 *(k <·> a2))) :+ 0.0+ phase01 _ _ = error "KagomeLattice: wrong dimensionality"+ phase12 k (a1:_:_) = (2 * cos(0.5 *(k <·> a1))) :+ 0.0+ phase12 _ _ = error "KagomeLattice: wrong dimensionality"+ phase20 k (a1:a2:_) = (2 * cos(0.5 *(k <·> (a1-a2)))) :+ 0.0+ phase20 _ _ = error "KagomeLattice: wrong dimensionality"++afOrder :: Parameterized Hamiltonian+afOrder = do j <- getScalar "J"+ d0 <- getVector "d0"+ d1 <- getVector "d1"+ d2 <- getVector "d2"+ return $ \_ -> scale (negate j)+ $ fromBlocks [[d0 .* s , 0 , 0 ]+ ,[0 , d1 .* s , 0 ]+ ,[0 , 0 , d2 .* s ]]+ where s = [sigmaX, sigmaY, sigmaZ]+ (.*) d ss = sum $ zipWith scale (toList d) ss++spinOrbit :: Parameterized Hamiltonian+spinOrbit = do t <- getScalar "tSO"+ as <- primitiveLattice+ return $ \k -> plusCT+ $ scale (iC * t)+ $ fromBlocks [[ 0 , scale (p01 k as) $ n01 .* s , 0 ]+ ,[ 0 , 0 , scale (p12 k as) $ n12 .* s ]+ ,[ scale (p20 k as) $ n20 .* s , 0 , 0 ]]++ where n01 = fromList [ cos ang01 , sin ang01 , 0.0 ]+ n12 = fromList [ cos ang12 , sin ang12 , 0.0 ]+ n20 = fromList [ cos ang20 , sin ang20 , 0.0 ]++ p01 k (_:a2:_) = (2 * cos(0.5 *(k <·> a2))) :+ 0.0+ p01 _ _ = error "KagomeLattice: wrong dimensionality"+ p12 k (a1:_:_) = (2 * cos(0.5 *(k <·> a1))) :+ 0.0+ p12 _ _ = error "KagomeLattice: wrong dimensionality"+ p20 k (a1:a2:_) = (2 * cos(0.5 *(k <·> (a1-a2)))) :+ 0.0+ p20 _ _ = error "KagomeLattice: wrong dimensionality"++ s = [sigmaX, sigmaY, sigmaZ]+ (.*) d ss = sum $ zipWith scale (toList d) ss+ plusCT m = m + tr m++ ang01 = pi/2.0 + 4.0*pi/3.0+ ang12 = pi/2.0+ ang20 = pi/2.0 + 2.0*pi/3.0
+ TBit/Systems/SquareLattice.hs view
@@ -0,0 +1,20 @@+module TBit.Systems.SquareLattice where++import TBit.Types+import TBit.Parameterization+import Numeric.LinearAlgebra.HMatrix+import Data.Complex (conjugate)+import Data.Map as M hiding ((!))++parameters = Parameters { scalarParams = loadParams [("m", 1.0)]+ , vectorParams = M.empty+ , latticeData = [ vector [0,1] , vector [1,0] ]+ , meshingData = Spacing 0.025}++hgteHamiltonian :: Parameterized Hamiltonian+hgteHamiltonian = do m <- getScalar "m"+ let mass k = 2 + m - ((cos (k!0) + cos(k!1)) :+ 0.0)+ let off k = sin (k!0) :+ (negate $ sin (k!1))+ return $ \k -> (2 >< 2)+ $ [ mass k , off k+ , conjugate (off k) , negate $ mass k ]
+ TBit/Toolbox.hs view
@@ -0,0 +1,56 @@+{-|+Module : TBit.Toolbox+Description : re-exports useful functions+Copyright : (c) Matthew Daniels, 2014+License : New BSD+Maintainer : danielsmw@gmail.com+Stability : experimental++This module re-exports several useful functions so that we don't need+to import a bunch of libraries at once into a script we'll frequently+be switching the output from.+-}++module TBit.Toolbox ( + -- * Topological calculations+ chern+ , chernBand+ , bandCurvature+ , occupiedCurvature+ -- * Orbital magnetization functions+ , orbMag+ , intrinsicOM+ , bandIntrinsicOM+ -- * Electric properties+ , nernstConductivity+ -- * Energetic properties+ , eigenstates+ , eigenenergies+ -- * Some basic Hamiltonian tools+ , squareLattice+ , hexLattice+ , kagomeLattice+ , kagomeRibbon+ , ring+ , toMatrix+ , neighborTerm+ , onsiteTerm+ , parityStaggeredTerm+ , localMoments+ , rashbaZ+ -- * Plotting utilities+ , bzPlot+ , paramPlot+ , bandPlot+ , kPath ) where++import TBit.Topological.Chern (chern, chernBand)+import TBit.Topological.Curvature (bandCurvature, occupiedCurvature)+import TBit.Magnetic.OrbitalMagnetization (orbMag, intrinsicOM, bandIntrinsicOM)+import TBit.Electronic.Conductance+import TBit.Plots (bzPlot, paramPlot, bandPlot)+import TBit.Framework (kPath)+import TBit.Hamiltonian.Eigenstates (eigenenergies, eigenstates)+import TBit.Hamiltonian.Builder.Examples (squareLattice, ring, hexLattice, kagomeRibbon, kagomeLattice)+import TBit.Hamiltonian.Builder.Matrification (toMatrix)+import TBit.Hamiltonian.Builder.Terms (neighborTerm, onsiteTerm, parityStaggeredTerm, localMoments, rashbaZ)
+ TBit/Topological/Chern.hs view
@@ -0,0 +1,110 @@+module TBit.Topological.Chern (chern, chernRaw, chernBand) where++import TBit.Types+import TBit.Framework+import TBit.Hamiltonian.Eigenstates+import TBit.Topological.Curvature++import Numeric.LinearAlgebra.HMatrix hiding ((!))++import Prelude hiding (lookup)+import Data.Map hiding (map)+import Data.Maybe (fromJust, isJust)+import qualified Data.Traversable as T++import Control.Monad+import Control.Monad.Except (throwError)+++{-|+ Calculate the Chern number of the nth band (indexed from 0) by+ integrating the Berry curvature over the Brillouin zone. The 'BandIndex'+ parameter is passed directly to 'TBit.Topological.Curvature.bandCurvature',+ and should use the same conventions for specifying the band.++ The output is appropriately normalized by 1/2π. The integration is carried+ out using the TanhSinh quadrature method via+ 'TBit.Numerical.Integration.integrate'.+-}+chernBand :: BandIndex -> Hamiltonian -> Parameterized Chern+chernBand n h = liftM ((*) (0.5/pi)) $ bzIntegral $ bandCurvature n h++chernRaw :: BandIndex -> Hamiltonian -> Parameterized Chern+chernRaw n h = do wg <- waveGrid h n+ liftM sum . liftM elems + . T.sequence + . mapWithKey (\gid _ -> curvature wg n gid)+ $ wg++{-|+ Calculate the Chern number of the first n occupied bands+ by using a grid of closed loops and calculating many Berry phases+ using the discretized formula. This function is /guaranteed/ to+ return an integer result by rounding the actual calculation. It tries+ to determine if the Chern number is undefined due to a degeneracy, and+ if it is then it throws an error via the 'Control.Monad.Except.ExceptT'+ monad transformer.+-}+chern :: BandIndex -> Hamiltonian -> Parameterized Chern+chern n h = liftM (fromIntegral . round) $ chernRaw n h++waveGrid :: Hamiltonian + -> BandIndex + -> Parameterized (Grid [Eigenstate])+waveGrid h b = do bz <- meshBZ+ T.mapM eigvs bz+ where eigvs = (safelyTake b . eigensystem h)++safelyTake :: Int -> Parameterized [(Energy, Eigenstate)] -> Parameterized [Eigenstate]+safelyTake n msys = do sys <- msys+ case compare n (length sys)+ of GT -> throwError dimErr+ EQ -> return . map snd $ sys+ LT -> if (fst $ sys !! (n-1)) == (fst $ sys !! n)+ then throwError chernErr+ else return . map snd . take n $ sys+ where dimErr = DimensionalityError ("Cannot compute the Chern number for more "+ ++"energy bands than the system supports.")+ chernErr = UndefinedError ("The Chern number is undefined due to a band "+ ++"degeneracy.")++curvature :: Grid [Eigenstate]+ -> BandIndex + -> GridIndex + -> Parameterized Double+curvature gr b (GID (m:n:[])) = do u1 <- phaseDiff gr b here right+ u2 <- phaseDiff gr b right there+ u3 <- phaseDiff gr b there up+ u4 <- phaseDiff gr b up here+ return $ ( phase+ $ u1 * u2 * u3 * u4 )+ / (2 * pi)+ where here = GID [m,n]+ right = GID [m+1,n]+ up = GID [m,n+1]+ there = GID [m+1,n+1]++curvature _ _ _ = throwError $ DimensionalityError + $ "Curvature calculations "+ ++ "are only implemented for 2D."++-- This will only work properly if there grid indices form a cube+phaseDiff :: Grid [Eigenstate]+ -> BandIndex+ -> GridIndex + -> GridIndex + -> Parameterized (Complex Double)+phaseDiff gr b m n = do let eigvs1 = lookup m gr+ let eigvs2 = lookup n gr++ let e1 = fromJust eigvs1+ let e2 = fromJust eigvs2++ case (isJust eigvs1, isJust eigvs2)+ of (True , True) -> return . det . (b >< b) . concat+ $ map (\v -> map (dot v) e2) e1+ (True , False) -> phaseDiff gr b ( m) (cyc n)+ (False, True) -> phaseDiff gr b (cyc m) ( n)+ (False, False) -> phaseDiff gr b (cyc m) (cyc n)++ where cyc (GID xs) = GID (map (\x -> if x == maximum xs then 0 else x) xs)
+ TBit/Topological/Curvature.hs view
@@ -0,0 +1,83 @@+module TBit.Topological.Curvature where ++import TBit.Types+import TBit.Framework+import TBit.Parameterization+import TBit.Numerical.Derivative+import TBit.Hamiltonian.Eigenstates++import Control.Monad++import Data.Map (elems)+import Data.List.Stream (delete)+import qualified Data.Traversable as T++import Numeric.LinearAlgebra.HMatrix++{-|+ Calculate the Berry curvature of a single band, which is to be given+ indexed from zero (i.e. to calculate the lowest band, pass in 0 for+ the 'BandIndex'. Uses the five-point stencil method for differentiation.+-}+bandCurvature :: BandIndex -> Hamiltonian -> Wavevector -> Parameterized Curvature+bandCurvature n h k = do ket <- eigenkets h k+ bra <- eigenbras h k+ eng <- eigenenergies h k+ return . negate . twice . imagPart . sum+ $ [ (num $ (bra!!n) <> hx <> (ket!!m))+ * (num $ (bra!!m) <> hy <> (ket!!n))+ / ((eng!!m - eng!!n)^2 :+ 0.0)+ | m <- delete n [0..pred dim]]+ where num m = m ! 0 ! 0+ dim = fst $ size $ h k+ twice = (*) 2.0+ kx = k ! 0+ ky = k ! 1+ hx = diff 0.0005 (\x -> h $ vector [realPart x, ky]) (kx :+ 0.0)+ hy = diff 0.0005 (\y -> h $ vector [kx, realPart y]) (ky :+ 0.0)++{-|+ Calculate the total Berry curvature of a the occupied bands, which are+ specified by passing in the number of filled bands as the first argument.+ For example, to find the curvature due to occupied bands of a 4 band system+ at half-filling, pass in 2 for the 'BandIndex'. Uses the five-point stencil+ method for differentiation.+-}+occupiedCurvature :: BandIndex -> Hamiltonian -> Wavevector -> Parameterized Curvature+occupiedCurvature b h k = do ket <- eigenkets h k+ bra <- eigenbras h k+ eng <- eigenenergies h k+ (Spacing s) <- getMesh+ return $ negate . twice . imagPart . sum+ $ zipWith (/)+ [ num $ (bra!!n) <> hx s <> (ket!!m)+ <> (bra!!m) <> hy s <> (ket!!n)+ | m <- occ , n <- unocc ]+ [ (eng!!m - eng!!n)^2 :+ 0.0+ | m <- occ , n <- unocc ]+ where num m = m ! 0 ! 0+ occ = [0 .. pred b]+ unocc = [b .. pred dim]+ dim = fst $ size $ h k+ twice = (*) 2.0+ kx = k ! 0+ ky = k ! 1+ hx s = diff (s :+ 0.0) (\x -> h $ vector [realPart x, ky]) (kx :+ 0.0)+ hy s = diff (s :+ 0.0) (\y -> h $ vector [kx, realPart y]) (ky :+ 0.0)+++-- |Deprecated?+curvatureFieldBand :: BandIndex -> Hamiltonian -> Parameterized [(Wavevector, Curvature)]+curvatureFieldBand n h = do grid <- meshBigBZ+ let field = fmap (\k -> do { bc <- bandCurvature n h k;+ return (k,bc) }) grid+ liftM elems $ T.sequence field+ +-- |Deprecated?+curvatureFieldOcc :: BandIndex -> Hamiltonian -> Parameterized [(Wavevector, Curvature)]+curvatureFieldOcc n h = do grid <- meshBZ+ let field = fmap (\k -> do { bc <- occupiedCurvature n h k;+ return (k,bc) }) grid+ liftM elems $ T.sequence field++
+ TBit/Topological/Spillage.hs view
@@ -0,0 +1,49 @@+{-|+Module : TBit.Topological.Spillage+Description : calculate spin-orbit and other spillages+Copyright : (c) Matthew Daniels, 2014+License : New BSD+Maintainer : danielsmw@gmail.com+Stability : experimental++Provide routines for determining the spillage, and in particular+Vanderbilt and Liu's spin-orbit spillage, given a Hamiltonian whose+spin-orbit coupling term can be turned on and off.+-}+module TBit.Topological.Spillage (spillage) where++import TBit.Types+import TBit.Parameterization (crunch, getScalar)+import TBit.Hamiltonian.Eigenstates (eigenbras,eigenkets)+import Debug.Trace++import Numeric.LinearAlgebra.HMatrix ((<>), fromBlocks, size, (!), Matrix)++import Control.Monad (liftM)+import Control.Monad.State++import Data.Map (adjust)+import Data.Complex++ +spillage :: String -> Filling -> Parameterized Hamiltonian -> Wavevector -> Parameterized Double+spillage ps n m k = do h' <- (modify (setScalar ps 0) >> m)+ h <- m+ s <- spillMatrix h h' k+ return $ (fromIntegral n) + - sum [ (realPart $ abs $ s!j!k)^2 + | j <- [0 .. pred n]+ , k <- [0 .. pred n]]+ where setScalar str val ps = ps { scalarParams = adjust (\_ -> val :+ 0.0) str + $ scalarParams ps }+ +spillMatrix :: Hamiltonian + -> Hamiltonian + -> Wavevector+ -> Parameterized (Matrix (Complex Double))+spillMatrix h h' k = do es <- eigenkets h k+ es' <- eigenbras h' k+ let dim = fst . size . h $ k+ return $ fromBlocks+ $ [[ (es'!!m) <> (es!!n) | m <- [0 .. pred dim]]+ | n <- [0 .. pred dim]]
+ TBit/Types.hs view
@@ -0,0 +1,165 @@+module TBit.Types ( Parameterized+ , Parameters (Parameters)+ , latticeData , meshingData+ , scalarParams , vectorParams+ , decomData+ , Lattice+ , TBError (..)+ , Grid+ , BandIndex+ , Filling+ , Meshing (Spacing)+ , GridIndex (GID)+ , Chern+ , Curvature+ , Wavevector+ , KPath+ , Hamiltonian+ , Magnetization+ , Moment+ , ChemEnergy+ , Energy+ , Eigenstate+ , Eigenbra+ , Eigenket+ , AFOrder + , Hopping + , OnSite + , Rashba + , SOC+ , Parameterizable+ , SiteData (..)+ , Displacement+ , CellGraph+ , AdjMatrix+ , Term+ , sigmaX+ , sigmaY+ , sigmaZ+ ) where++import Control.Applicative+import Control.DeepSeq+import Data.Functor()+import Data.Foldable+import Data.Traversable+import Data.Map (Map)+import Control.Monad.Except (ExceptT)+import Control.Monad.State (State)+import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, iC, (><))+import Data.Complex (Complex (..))+import Data.Graph.Inductive hiding ((><))+import Control.DeepSeq ()++-- Computational types+type Parameterized = ExceptT TBError (State Parameters)+data Parameters = Parameters { latticeData :: Lattice+ , meshingData :: Meshing+ , decomData :: [LEdge Displacement]+ , scalarParams :: Map String (Complex Double)+ , vectorParams :: Map String (Vector (Complex Double)) }+ deriving (Show)++type Term = String -> CellGraph -> Parameterized (Wavevector -> AdjMatrix)++instance NFData Meshing where+ rnf (Spacing t) = t `seq` ()++instance NFData Parameters where+ rnf ps = map rnf (latticeData ps) + `seq` rnf (meshingData ps)+ `seq` map rnf (decomData ps)+ `seq` rnf (scalarParams ps)+ `seq` rnf (vectorParams ps)++type Lattice = [Vector Double]+data TBError = SingularLatticeError + | DimensionalityError String+ | UndefinedError String+ | UnknownParameter String+ deriving Show++type Grid = Map GridIndex+type BandIndex = Int+type Filling = BandIndex++data Meshing = Spacing Double deriving Show+data GridIndex = GID [Int] deriving (Ord, Eq, Show)++-- Physical quantities+type Chern = Double+type Wavevector = Vector Double+type Curvature = Double+type KPath = [(String,Wavevector)]+type Hamiltonian = Wavevector -> Matrix (Complex Double)+type Magnetization = Double+type Moment = Vector Double++type ChemEnergy = Double+type Energy = Double+type Eigenstate = Vector (Complex Double)+type Eigenbra = Matrix (Complex Double)+type Eigenket = Matrix (Complex Double)++-- TB Parameters+type AFOrder = Complex Double+type Hopping = Complex Double+type OnSite = Complex Double+type Rashba = Complex Double+type SOC = Complex Double+type Position = Vector Double+data SiteData = ScalarSite { num :: Int }+ | VectorSite { num :: Int, mom :: Moment } deriving Show++type Parameterizable = ExceptT TBError (State Parameters)+type Displacement = Vector Double+type CellGraph = Gr SiteData Displacement+type AdjMatrix = Gr (Matrix (Complex Double)) (Matrix (Complex Double))++instance Traversable ((,) a) where+ traverse f (x,z) = (,) x <$> f z++instance Traversable ((,,) a b) where+ traverse f (x,y,z) = (,,) x y <$> f z++instance Foldable ((,) a) where+ foldMap f (_, z) = f z+ foldr f t (_, z) = f z t++instance Foldable ((,,) a b) where+ foldMap f (_, _, z) = f z+ foldr f t (_, _, z) = f z t++instance Functor ((,,) a b) where+ fmap f (x,y,z) = (x,y,f z)++instance Traversable ((,,,) a b c) where+ traverse f (x,y,z,w) = (,,,) x y z <$> f w++instance Foldable ((,,,) a b c) where+ foldMap f (_, _, _, w) = f w+ foldr f t (_, _, _, w) = f w t++instance Functor ((,,,) a b c) where+ fmap f (x,y,z,w) = (x,y,z,f w)++instance NFData TBError where+ rnf e@(SingularLatticeError) = e `seq` ()+ rnf (DimensionalityError s) = s `seq` ()+ rnf (UndefinedError s) = s `seq` ()+ rnf (UnknownParameter s) = s `seq` ()++sigmaX :: Matrix (Complex Double)+sigmaX = (2 >< 2)+ [ 0 , 1+ , 1 , 0 ]++sigmaY :: Matrix (Complex Double)+sigmaY = (2 >< 2)+ [ 0 , -iC+ , iC , 0 ]++sigmaZ :: Matrix (Complex Double)+sigmaZ = (2 >< 2)+ [ 1 , 0+ , 0 ,-1 ]