RFIM.hs

import System.Environment (getArgs)
import System.Random
import Control.Monad      (liftM)
import Data.Array.Unboxed
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU

configs :: Int -> [[Char]] -- List of configurations
configs n = V.toList $ V.take (3^n) $ V.fromList $ iterate nextConfig $ replicate n '0'
        where nextConfig ('0':xs) = '1':xs
              nextConfig ('1':xs) = '2':xs
              nextConfig ('2':xs) = '0':(nextConfig xs)

normals :: [Double] -> [Double] -- I.i.d. N(0,1) random variables
normals (r1:r2:rs) = (sqrt (-2 * log r1) * cos (2*pi * r2)) : (normals rs)

data Stats = Stats { n0 :: !Int, n1 :: !Int, n2 :: !Int, ff :: !Double }

curieWeiss :: [Double] -> [Char] -> (Int,Double) -- Field -> Config -> RF Hamiltonian
curieWeiss f xs = (n1-n2, exp(field-det))
           where Stats n0 n1 n2 ff = VU.foldl acc (Stats 0 0 0 0) $ VU.zip (VU.fromList xs) (VU.fromList f)
                 acc s@(Stats k0 _ _ uu) ('0',u) = s { n0 = k0 + 1, ff = uu +     u }
                 acc s@(Stats _ k1 _ uu) ('1',u) = s { n1 = k1 + 1, ff = uu - 0.5*u }
                 acc s@(Stats _ _ k2 uu) ('2',u) = s { n2 = k2 + 1, ff = uu - 0.5*u }
                 field = (fromIntegral $ n0+n1+n2) * ff
                 det   = fromIntegral $ 3 * (n0*n1 + n0*n2 + n1*n2)

run :: [Double] -> Int -> UArray Int Double -- Field -> n -> partition function, fractioned by n1-n2
run f n = accumArray (+) 0.0 (-n,n) $ map (curieWeiss f) $ configs n

main :: IO ()
main = do n <- liftM (read . head) getArgs
          f <- liftM (take n . normals . randoms) getStdGen

          mapM_ p $ assocs $ run (map sign f) n
                where p (i,x) = putStrLn $ (show i) ++ " " ++ (show $ - log x)
                      sign x  = if x>0 then 1.0 else -1.0