Adding the split-operator chapter and code (algorithm-archivists#108)

* adding smallscale changes.

* adding partial draft for the split-op method

* adding draft of split-op chapter.

* adding back split-op.jl file.

* fixing typo in split-op method.

* turning the splitop param module into a struct

* adding smallscale changes to spacing in split_op.jl

Author: leios

README.md

@@ -4,7 +4,9 @@ This goal is obviously too ambitious for a book of any size, but it is a great p
 The book can be found here: https://www.algorithm-archive.org/.
 The github repository can be found here: https://github.com/algorithm-archivists/algorithm-archive.
 Most algorithms have been covered on the youtube channel LeiosOS: https://www.youtube.com/user/LeiosOS
-and livecoded on Twitch: https://www.twitch.tv/simuleios
+and livecoded on Twitch: https://www.twitch.tv/simuleios.
+If you would like to communicate more directly, please feel free to go to our discord: https://discord.gg/Pr2E9S6.
+
 
 Note that the this project is is essentially a book about algorithms collaboratively written by an online community. 
 Fortunately, there are a lot of algorithms out there, which means that there is a lot of content material available.

chapters/physics_solvers/quantum/quantum.md

@@ -57,6 +57,16 @@ $$
 \int_{-\infty}^{+\infty}|\Psi(\mathbf{r},t)|^2 d\mathbf{r} = 1
 $$
 
+As another note: Just like position space can be parameterized by a position vector $$\textbf{x}$$, wavefunctions can be parameterized by a _wave_vector $$\textbf{k}$$ in frequency space.
+Often times, the wavevector space is called _momentum_ space, which makes sense when considering the de Broglie formula:
+
+$$
+p = \frac{h}{\lambda}
+$$
+
+where $$h$$ is Planck's constant and $$\lambda$$ is the wavelength.
+This means that we can ultimately move between real and momentum space by using [Fourier Transforms](../../FFT/cooley_tukey.md), which is incredibly useful in a number of cases!
+
 Unfortunately, the interpretation of quantum simulation is rather tricky and is ultimately easier to understand with slightly different notation. 
 This notation is called _braket_ notation, where a _ket_ looks like this:
 

chapters/physics_solvers/quantum/split-op/code/julia/split_op.jl

@@ -0,0 +1,79 @@
+using Plots
+pyplot()
+
+# struct to hold all parameters for simulation
+struct Param
+    xmax::Float64
+    res::Int64
+    dt::Float64
+    timesteps::Int64
+    dx::Float64
+    x::Vector{Float64}
+    dk::Float64
+    k::Vector{Float64}
+
+    Param() = new(10.0, 512, 0.05, 1000, 2 * 10.0/512,
+                  Vector{Float64}(-10.0 + 10.0/512 : 20.0/512 : 10.0),
+                   pi / 10.0,
+                  Vector{Float64}(vcat(0:512/2 - 1, -512/2 : -1) * pi/10.0))
+    Param(xmax::Float64, res::Int64, dt::Float64, timesteps::Int64) = new(
+              xmax, res, dt, timesteps, 
+              2*xmax/res, Vector{Float64}(-xmax+xmax/res:2*xmax/res:xmax),
+              pi/xmax, Vector{Float64}(vcat(0:res/2-1, -res/2:-1)*pi/xmax)
+          )
+end
+
+# struct to hold all operators
+mutable struct Operators
+    V::Vector{Complex{Float64}}
+    PE::Vector{Complex{Float64}}
+    KE::Vector{Complex{Float64}}
+    wfc::Vector{Complex{Float64}}
+end
+
+# Function to initialize the wfc and potential
+function init(par::Param, voffset::Float64, wfcoffset::Float64)
+    V = 0.5 * (par.x - voffset).^2
+    wfc = 3 * exp.(-(par.x - wfcoffset).^2/2)
+    PE = exp.(-0.5*im*V*par.dt)
+    KE = exp.(-0.5*im*par.k.^2*par.dt)
+
+    opr = Operators(V, PE, KE, wfc)
+end
+
+# Function for the split-operator loop
+function split_op(par::Param, opr::Operators)
+
+    for i = 1:par.timesteps
+        # Half-step in real space
+        opr.wfc = opr.wfc .* opr.PE
+
+        # fft to phase space
+        opr.wfc = fft(opr.wfc)
+
+        # Full step in phase space
+        opr.wfc = opr.wfc .* opr.KE
+
+        # ifft back
+        opr.wfc = ifft(opr.wfc)
+
+        # final half-step in real space
+        opr.wfc = opr.wfc .* opr.PE
+
+        # plotting density and potential
+        density = abs2.(opr.wfc)
+
+        plot([density, real(opr.V)])
+        savefig("density" * string(lpad(i, 5, 0)) * ".png")
+        println(i)
+    end
+end
+
+# main function
+function main()
+    par = Param(10.0, 512, 0.05, 1000)
+    opr = init(par, 0.0, 1.0)
+    split_op(par, opr)
+end
+
+main()

chapters/physics_solvers/quantum/split-op/split-op.md

@@ -17,9 +17,76 @@ This is the system I studied for most of my PhD (granted, we played a few tricks
 
 At it's heart, the split-op method is nothing more than a pseudo-spectral differential equation solver... That is to say, it solves the Schrodinger equation with [FFT's](../../../FFT/cooley_tukey.md).
 In fact, there is a large class of spectral and pseudo-spectral methods used to solve a number of different physical systems, and we'll definitely be covering those in the future.
+As mentioned in the [quantum systems](../quantum.md) section, we can represent a a quantum wavefunction in momentum space, which is parameterized with the wavevector $$k$$.
+In the hamiltonian shown above, we can split our system into real-space components, $$\hat{H}_R = \left[V(\mathbf{r}) + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t)$$, and momentum space components, $$\hat{H}_M = \left[-\frac{\hbar^2}{2m}\nabla^2 \right]\Psi(\mathbf{r},t)$$.
+If we assume a somewhat general solution to our quantum system:
 
-For now, let's answer the obvious question, "How on Earth can FFT's be used to solve quantum systems?"
-That is precisely the question we will answer here!
+$$
+\Psi(\mathbf{r},t + dt) = \left[e^{-\frac{i\hat{H}dt}{\hbar}}\right]\Psi(\mathbf{r},t) = \left[e^{-\frac{i(\hat{H}_R + \hat{H}_M)dt}{\hbar}}\right]\Psi(\mathbf{r},t)
+$$
+
+and assume we are simulating our system by a series of small tiemsteps ($$dt$$), we can perform similar splitting by using the Baker-Campbell-Housdorff formula:
+
+$$
+\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_Rdt}{\hbar}}e^{-\frac{i\hat{H}_Mdt}{\hbar}}e^{-\frac{[i\hat{H}_R, i\hat{H}_M]dt^2}{2}}\right]\Psi(\mathbf{r},t)
+$$
+
+This accrues a small amount of error ($$dt^2$$) related to the commutation of the real and momentum-space components of the Hamiltonian. That's not okay.
+In order to change the $$dt^2$$ error to $$dt^3$$, we can split the system by performing a half-step in real space before doing a full-step in momentum space, through a process called _Strang Splitting_ like so:
+
+$$
+\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_Rdt}{2\hbar}}e^{-\frac{i\hat{H}_Mdt}{\hbar}}e^{-\frac{i\hat{H}_Rdt}{2\hbar}} \right]\Psi(\mathbf{r},t) + \mathcal{O}(dt^3)
+$$
+
+We can then address each part of this solution in chunks, first in real space, then in momentum space, then in real space again by using [Fourier Transforms](../../../FFT/cooley_tukey.md).
+Which looks something like this:
+
+$$
+\Psi(\mathcal{r}, t+dt) = \left[\hat{U}_R(\frac{dt}{2})\mathcal{F}\left[\hat{U}_M(dt) \mathcal{F} \left[\hat{U}_R(\frac{dt}{2}) \Psi(\mathbf{r},t) \right] \right] \right] + \mathcal{O}(dt^3)
+$$
+
+where $$\hat{U}_R = e^{-\frac{i\hat{H}_Rdt}{\hbar}}$$, $$\hat{U}_M = e^{-\frac{i\hat{H}_Mdt}{\hbar}}$$, and $$\mathcal{F}$$ and $$\mathcal{F}^{-1}$$ indicate forward and inverse Fourier Transforms.
+
+As a small concession here, using this method enforces periodic boundary conditions, where the wavefunction will simply slide from one side of your simulation box to the other, but that's fine for most cases.
+In fact, for many cases (such as large-scale turbulence models) it's ideal.
+
+Luckily, the code in this case is pretty straightforward.
+Frist, we need to set all the initial parameters, including the initial grids in real and momentum space:
+
+{% method %}
+{% sample lang="jl" %}
+[import:4-31, lang:"julia"](code/julia/split_op.jl)
+{% 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`.
+This is simply because the FFT will naturally assume that the `0` in our grid is at the left side of the simulation, so we shift k-space to match this expectation.
+Afterwards, we turn them into operators:
+
+{% method %}
+{% sample lang="jl" %}
+[import:32-41, lang:"julia"](code/julia/split_op.jl)
+{% 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.
+As a note, if we run this simulation in _imaginary time_, by simply setting $$\tau = it$$ and stepping through $$\tau$$, we will no longer see an "real-world" example of how the atoms should behave, but will instead see an exponential decay of higher-energy states.
+This means that we can find the ground state of our system by running the simulation in imaginary time, which is an incredibly useful feature!
+
+And finally go step-by-step through the simulation:
+
+{% method %}
+{% sample lang="jl" %}
+[import:42-69, lang:"julia"](code/julia/split_op.jl)
+{% endmethod %}
+
+And that's it.
+The Split-Operator method is one of the most commonly used quantum simulation algorithms because of how straightforward it is to code and how quickly you can start really digging into the physics of the simulation results!
+
+# Example Code
+{% method %}
+{% sample lang="jl" %}
+### Julia
+[import, lang:"julia"](code/julia/split_op.jl)
+{% endmethod %}
 
 <script>
 MathJax.Hub.Queue(["Typeset",MathJax.Hub]);

update_site.sh

@@ -26,7 +26,7 @@ if [ $# -ge 1 ]; then
     cd "$out_dir"
     git reset --hard && \
       git clean -dfx && \
-      git pull
+      git pull origin master
 
     if [ $? -ne 0 ]; then
         echo "Failed to prepare repository"