Split-op in Haskell (#357)

* Added split-op in Haskell.
Added energy calculations.
Added a gnuplotscript to visualize wavefunctions.

* Changed contributors? Not sure

* Included code to chapter

* Update all julia code to V1.0 (#389)

* adding changes to update all julia code to V1.0

* adding more points to graham.jl

* Verlet integration in Fortran90 (#364)

* initial commit

* Programming done

* Minor style changes. Merge ready

* hopefully fix PR

* minor edits for code readability

* indexed main procedure

* fix typo

* Euclidean's algorithm in Fortran90 (#369)

* implementation and documentary

* minor code clean up

* fixing git again...

* Delete monte_carlo_pi.f90

again fortan file automagically appears.

* code style changes

* code style

* indexed main procedure

* Cleaned up code

* Refactored .*, .+ and normalization operators
* Got rid of {-# LANGUAGE FlexibleContexts #-} by giving explicit signature to normalize
* Small cosmetic changes

* Undid change to CONTTRIBUTORS.md

* Moved gnuplot script to own folder and updated instructions to use it

* Removed old plot script

* Changed import lines

Author: jiegillet

contents/quantum_systems/code/haskell/Energy.hs

@@ -0,0 +1,14 @@
+import Data.Array.CArray
+import Data.Complex
+import Math.FFT (dft, idft) -- Binding to fftw
+
+type Vector = CArray Int (Complex Double)
+
+calculateEnergy :: Double -> Vector -> Vector -> Vector -> Double
+calculateEnergy dx kin pot wfc = (* dx) . sum . map realPart $ elems total
+  where
+    total = liftArray2 (+) kineticE potentialE
+    potentialE = wfcConj .* pot .* wfc
+    kineticE = wfcConj .* idft (kin .* dft wfc)
+    wfcConj = liftArray conjugate wfc
+    a .* b = liftArray2 (*) a b

contents/quantum_systems/quantum_systems.md

@@ -227,6 +227,8 @@ This ultimately looks like this:
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/energy.jl)
+{% sample lang="hs" %}
+[import, lang:"haskell"](code/haskell/Energy.hs)
 {% sample lang="c" %}
 [import:28-, lang:"c_cpp"](code/c/energy.c)
 {% sample lang="py" %}

contents/split-operator_method/code/gnuplot/plot_output.plt

@@ -0,0 +1,24 @@
+# use like: gnuplot -e "folder='/path/to/data/directory'" plot_output.plt
+
+set terminal gif animate delay 10
+set output folder.'/wavefunction.gif'
+
+set key off
+
+set xrange [-10:10] 
+set yrange [0:1]
+set y2range [0:50]
+
+set ytics nomirror autofreq tc lt 1
+set ylabel 'Psi(x)' tc lt 1
+
+set y2tics nomirror autofreq tc lt 2
+set y2label 'V(x)' tc lt 2
+
+list = system('ls '.folder.'/output*.dat')
+
+do for [file in list] {
+    plot file u 1:2 smooth csplines, \
+         file u 1:3 smooth csplines axes x1y2
+}
+

contents/split-operator_method/code/haskell/splitOp.hs

@@ -0,0 +1,104 @@
+import Data.Array.CArray
+import Data.Complex
+import Data.List (intercalate, transpose)
+import Math.FFT (dft, idft)
+
+type Vector = CArray Int (Complex Double)
+
+(.*), (.+) :: Vector -> Vector -> Vector
+a .* b = liftArray2 (*) a b
+a .+ b = liftArray2 (+) a b
+
+normalize :: Double -> Vector -> Vector
+normalize dx v =
+  let factor = 1 / sqrt dx / norm2 v :+ 0
+   in liftArray (factor *) v
+
+data Parameters = Parameters
+  { xmax :: Double
+  , res :: Int
+  , dt :: Double
+  , timesteps :: Int
+  , dx :: Double
+  , x :: Vector
+  , dk :: Double
+  , ks :: Vector
+  , imTime :: Bool
+  }
+
+defaultParameters :: Parameters
+defaultParameters = makeParameters 10 512 0.01 1000 True
+
+makeParameters :: Double -> Int -> Double -> Int -> Bool -> Parameters
+makeParameters xmax res dt timesteps imTime =
+  let fi = fromIntegral
+      rng = (0, res - 1)
+      ks = [0 .. div res 2 - 1] ++ [-div res 2 .. -1]
+   in Parameters
+        xmax
+        res
+        dt
+        timesteps
+        (2 * xmax / fi res)
+        (listArray rng $
+         map (\n -> xmax * (-1 + 2 * fi n / fi res) :+ 0) [1 .. res])
+        (pi / xmax)
+        (listArray rng $ map ((:+ 0) . (pi / xmax *) . fi) ks)
+        imTime
+
+data Operators = Operators
+  { v :: Vector
+  , rStep :: Vector
+  , kStep :: Vector
+  , wfc :: Vector
+  }
+
+makeOperators :: Parameters -> Complex Double -> Complex Double -> Operators
+makeOperators param v0 wfc0 =
+  let rng = (0, res param - 1)
+      time
+        | imTime param = dt param :+ 0
+        | otherwise = 0 :+ dt param
+      v = liftArray (\x -> 0.5 * (x - v0) ^ 2) (x param)
+      rStep = liftArray (\x -> exp (-0.5 * time * x)) v
+      kStep = liftArray (\k -> exp (-0.5 * time * k ^ 2)) (ks param)
+      wfc = liftArray (\x -> exp (-(x - wfc0) ^ 2 / 2)) (x param)
+   in Operators v rStep kStep (normalize (dx param) wfc)
+
+evolve :: Parameters -> Operators -> [Operators]
+evolve param op@(Operators _ rStep kStep _) = iterate splitop op
+  where
+    splitop op = op {wfc = wfc' op}
+    wfc' = norm . (rStep .*) . idft . (kStep .*) . dft . (rStep .*) . wfc
+    norm = if imTime param then normalize (dx param) else id
+
+calculateEnergy :: Parameters -> Operators -> Double
+calculateEnergy param ops = (* dx param) . sum . map realPart $ elems totalE
+  where
+    totalE = potentialE .+ kineticE
+    potentialE = wfcConj .* v ops .* wfc ops
+    kineticOp = liftArray ((/ 2) . (^ 2)) (ks param)
+    kineticE = wfcConj .* idft (kineticOp .* dft (wfc ops))
+    wfcConj = liftArray conjugate $ wfc ops
+
+-- Use gnuplot to make an animated  GIF using ../gnuplot/plot_output.plt
+-- $ gnuplot -e "folder='../haskell'" plot_output.plt
+printEvolution :: Parameters -> [Operators] -> IO ()
+printEvolution param =
+  mapM_ (export . (format <$>)) . zip [0 ..] . take 100 . skip
+  where
+    skip (x:xs) = x : skip (drop (div (timesteps param) 100 - 1) xs)
+    format (Operators v _ _ wfc) =
+      let density = liftArray ((^ 2) . abs) wfc
+          values = map (map (show . realPart) . elems) [x param, density, v]
+       in intercalate "\n" $ map (intercalate "\t") $ transpose values
+    export (i, f) = writeFile ("output" ++ pad (show i) ++ ".dat") f
+    pad n = replicate (5 - length n) '0' ++ n
+
+main :: IO ()
+main = do
+  let p = defaultParameters
+      o = makeOperators p 0 4
+      evol = evolve p o
+  print $ calculateEnergy p (evol !! timesteps p)
+  printEvolution p evol

contents/split-operator_method/split-operator_method.md

@@ -55,7 +55,7 @@ Here's a flowchart of what we are looking for every timestep:
 
 
 For the most part, that's it:
-1. Multiply the wavefunction in real space with the real-space operator. 
+1. Multiply the wavefunction in real space with the real-space operator.
 2. Flip to momentum space with a Fourier transform.
 3. Multiply the momentum-space wavefuntion by the momentum-space operator.
 4. Flip to position space with an inverse Fourier transform.
@@ -104,6 +104,8 @@ Regardless, we first need to set all the initial parameters, including the initi
 [import:51-72, lang:"c_cpp"](code/c/split_op.c)
 {% sample lang="py" %}
 [import:11-30, lang:"python"](code/python/split_op.py)
+{% sample lang="hs" %}
+[import:17-47, lang:"haskell"](code/haskell/splitOp.hs)
 {% endmethod %}
 
 As a note, when we generate our grid in momentum space `k`, we need to split the grid into two lines, one that is going from `0` to `-kmax` and is then discontinuous and goes from `kmax` to `0`.
@@ -121,6 +123,8 @@ Afterwards, we turn them into operators:
 [import:74-95, lang:"c_cpp"](code/c/split_op.c)
 {% sample lang="py" %}
 [import:33-54, lang:"python"](code/python/split_op.py)
+{% sample lang="hs" %}
+[import:49-66, lang:"haskell"](code/haskell/splitOp.hs)
 {% endmethod %}
 
 Here, we use a standard harmonic potential for the atoms to sit in and a gaussian distribution for an initial guess for the probability distribution.
@@ -138,6 +142,8 @@ The final step is to do the iteration, itself.
 [import:97-145, lang:"c_cpp"](code/c/split_op.c)
 {% sample lang="py" %}
 [import:57-95, lang:"python"](code/python/split_op.py)
+{% sample lang="hs" %}
+[import:68-73, lang:"haskell"](code/haskell/splitOp.hs)
 {% endmethod %}
 
 And that's it.
@@ -161,6 +167,8 @@ Checking to make sure your code can output the correct energy for a harmonic tra
 [import, lang:"c_cpp"](code/c/split_op.c)
 {% sample lang="py" %}
 [import:5-127, lang:"python"](code/python/split_op.py)
+{% sample lang="hs" %}
+[import, lang:"haskell"](code/haskell/splitOp.hs)
 {% endmethod %}
 
 <script>