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