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