Prototype of Verlet method in Haskell

verlet.hs

import Data.List
import Control.Monad.Random
import Text.Printf
import System.Environment
import Debug.Trace

data PhaseSpace = PhaseSpace {
  particles :: [(Pos,Momentum)]
}

type D3 = (Double,Double,Double)
type Pos = D3
type Momentum = D3
type Force = D3

instance Show PhaseSpace where
  show (PhaseSpace ps) = "[" ++ (intercalate ", " $ map f ps) ++ "]"
    where 
      f ((rx,ry,rz),(px,py,pz)) = printf "r[%.3f %.3f %.3f]  p[%.3f %.3f %.3f]" rx ry rz px py pz

mytrace a = traceShow a a

rand3d :: (RandomGen g) => Double -> Rand g D3
rand3d a = do
  x <- getRandomR (-a,a)
  y <- getRandomR (-a,a)
  z <- getRandomR (-a,a)
  return (x,y,z)

getInitial :: Int -> IO PhaseSpace
getInitial n = do
  ps <- evalRandIO (ps_seq 1000)
  return $ PhaseSpace $ zip rs ps
  where
    range = map fromIntegral [0..n-1]
    rs = [(x,y,z)| x <- range, y <- range, z <- range]
    ps_seq n = sequence (replicate n (rand3d 1))


-- This is slower than (filter pred . map f), since fusion transformation does not occur.
mapfilter :: (a -> b) -> (b -> Bool) -> [a] -> [b]
mapfilter _ _ [] = []
mapfilter f pred (x:xs)
  = let y = f x in
      if pred y then
        y:mapfilter f pred xs
      else
        mapfilter f pred xs

-- Van der Waals force
force :: [Pos] -> Int-> D3
force ps i = vsum $ map f $ g $ map ((ps!!i) -:) ps
  where
    --g :: [a] -> [a]
    g = filter ((> 0.1) . vnorm)
    g _ = filter (>0.1 . vnorm + 1)
    -- g rs = let (h, t) = splitAt i rs in h ++ (tail t)
    f (x,y,z) = (ff x y z,ff y z x,ff z x y)
    ff x y z = (a*2*x) * (r2**(-7)*(-6) + r2**(-4)*3)
      where r2 = x**2+y**2+z**2
    a = 0.01    

type Iso a = a -> a

infixl 3 `vplus`
infixl 3 +:
(+:) = vplus
vplus :: D3 -> D3 -> D3
vplus (a,b,c) (d,e,f) = (a+d,b+e,c+f)

infixl 3 -:
(a,b,c) -: (d,e,f) = (a-d,b-e,c-f)

vsum :: [D3] -> D3
vsum vs | null vs = (0,0,0)
        | otherwise = foldl1' (+:) vs

vnorm (x,y,z) = x**2 + y**2 + z**2

infixr 4 `vmul`
infixr 4 *:
(*:) = vmul
vmul :: Double -> D3 -> D3
vmul k (a,b,c) = (k*a,k*b,k*c)

type MoveFunc = Double -> Iso PhaseSpace
move :: MoveFunc
move dt (PhaseSpace particles) = PhaseSpace (zip rs_new ps_new)
  where
    ps = map snd particles
    ps_new :: [D3]
    ps_new = zipWith h [0..] ps
      where
        h :: Int -> Momentum -> D3
        h i p = p +: (dt/2) *: ((fs rs !! i) +: (fs rs_new !! i))
    rs :: [D3]
    rs = map fst particles
    rs_new :: [D3]
    rs_new = zipWith g [0..] rs
      where
        g :: Int -> Pos -> D3
        g i r = (r +: dt *: (ps!!i) +: dt**dt *: (fs rs !! i))
    fs :: [D3] -> [Force]
    fs rs_ = map (force rs_) [0..length rs_-1]

iterate' :: (a -> a) -> a -> [a]
iterate' f x = x `seq` (x : iterate' f (f x))

allForces (PhaseSpace particles)
  = map (force (map fst particles)) [0..length particles-1]

main = do
  as <- getArgs
  initial <- getInitial (read $ as !! 0)
  -- print $ allForces initial
  let final = (iterate' (move (read $ as !! 1)) initial) !! (read $ as !! 2)
  print final