Skip to content

Instantly share code, notes, and snippets.

@nebuta
Last active December 14, 2015 17:29
Show Gist options
  • Save nebuta/5122571 to your computer and use it in GitHub Desktop.
Save nebuta/5122571 to your computer and use it in GitHub Desktop.
Prototype of Verlet method in Haskell
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment