Split op updates

contents/quantum_systems/code/julia/energy.jl

@@ -0,0 +1,18 @@
+# We are calculating the energy to check <Psi|H|Psi>
+function calculate_energy(wfc, H_k, H_r, dx)
+    # Creating momentum and conjugate wavefunctions
+    wfc_k = fft(wfc_r)
+    wfc_c = conj(wfc_r)
+
+    # Finding the momentum and real-space energy terms
+    energy_k = wfc_c.*ifft((H_k) .* wfc_k)
+    energy_r = wfc_c.* H_r .* wfc_r
+
+    # Integrating over all space
+    energy_final = 0
+    for i = 1:length(energy_k)
+        energy_final += real(energy_k[i] + energy_r[i])
+    end 
+
+    return energy_final*dx
+end

contents/split-operator_method/code/julia/split_op.jl

@@ -1,7 +1,11 @@
-using Plots
-pyplot()
+#------------split_op.jl-------------------------------------------------------#
+#
+# Plotting: to plot individual timesteps, use gnuplot like so:
+#               p "output00000.dat" u 1:2 w l
+#               rep "output00000.dat" u 1:3 w l
+#
+#------------------------------------------------------------------------------#
 
-# struct to hold all parameters for simulation
 struct Param
     xmax::Float64
     res::Int64
@@ -11,34 +15,48 @@ struct Param
     x::Vector{Float64}
     dk::Float64
     k::Vector{Float64}
+    im_time::Bool
 
     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(
+                  pi / 10.0,
+                  Vector{Float64}(vcat(0:512/2 - 1, -512/2 : -1) * pi/10.0),
+                  false)
+    Param(xmax::Float64, res::Int64, dt::Float64, timesteps::Int64,
+          im_val::Bool) = 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)
+              pi/xmax, Vector{Float64}(vcat(0:res/2-1, -res/2:-1)*pi/(xmax)),
+              im_val
           )
 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}}
+
+    Operators(res) = new(Vector{Complex{Float64}}(res),
+                         Vector{Complex{Float64}}(res),
+                         Vector{Complex{Float64}}(res),
+                         Vector{Complex{Float64}}(res))
 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(length(par.x))
+    opr.V = 0.5 * (par.x - voffset).^2
+    opr.wfc = exp.(-(par.x - wfcoffset).^2/2)
+    if (par.im_time)
+        opr. = exp.(-0.5*par.k.^2*par.dt)
+        opr.PE = exp.(-0.5*opr.V*par.dt)
+    else
+        opr. = exp.(-im*0.5*par.k.^2*par.dt)
+        opr.PE = exp.(-im*0.5*opr.V*par.dt)
+    end
 
-    opr = Operators(V, PE, KE, wfc)
+    return opr
 end
 
 # Function for the split-operator loop
@@ -48,32 +66,78 @@ function split_op(par::Param, opr::Operators)
         # Half-step in real space
         opr.wfc = opr.wfc .* opr.PE
 
-        # fft to phase space
+        # fft to momentum space
         opr.wfc = fft(opr.wfc)
 
-        # Full step in phase space
-        opr.wfc = opr.wfc .* opr.KE
+        # Full step in momentum space
+        opr.wfc = opr.wfc .* opr.
 
         # ifft back
         opr.wfc = ifft(opr.wfc)
 
         # final half-step in real space
         opr.wfc = opr.wfc .* opr.PE
 
-        # plotting density and potential
+        # density for plotting and potential
         density = abs2.(opr.wfc)
 
-        plot([density, real(opr.V)])
-        savefig("density" * string(lpad(i, 5, 0)) * ".png")
-        println(i)
+        # renormalizing for imaginary time
+        if (par.im_time)
+            renorm_factor = sum(density) * par.dx
+
+            for j = 1:length(opr.wfc)
+                opr.wfc[j] /= sqrt(renorm_factor)
+            end
+        end
+
+        # Outputting data to file. Plotting can also be done in a similar way
+        # This is set to output exactly 100 files, no matter how many timesteps
+        if ((i-1) % div(par.timesteps, 100) == 0)
+            outfile = open("output" *string(lpad(i-1, 5, 0))* ".dat","w")
+
+            # Outputting for gnuplot. Any plotter will do.
+            for j = 1:length(density)
+                write(outfile, string(par.x[j]) * "\t"
+                               * string(density[j]) * "\t"
+                               * string(real(opr.V[j])) * "\n")
+            end
+
+            close(outfile)
+            println("Outputting step: ", i)
+        end
     end
 end
 
+# We are calculating the energy to check <Psi|H|Psi>
+function calculate_energy(par, opr)
+    # Creating real, momentum, and conjugate wavefunctions
+    wfc_r = opr.wfc
+    wfc_k = fft(wfc_r)
+    wfc_c = conj(wfc_r)
+
+    # Finding the momentum and real-space energy terms
+    energy_k = 0.5*wfc_c.*ifft((par.k.^2) .* wfc_k)
+    energy_r = wfc_c.*opr.V .* wfc_r
+
+    # Integrating over all space
+    energy_final = 0
+    for i = 1:length(energy_k)
+        energy_final += real(energy_k[i] + energy_r[i])
+    end 
+
+    return energy_final*par.dx
+end
+
 # main function
 function main()
-    par = Param(10.0, 512, 0.05, 1000)
-    opr = init(par, 0.0, 1.0)
+    par = Param(5.0, 256, 0.05, 100, true)
+
+    # Starting wavefunction slightly offset so we can see it change
+    opr = init(par, 0.0, -1.00)
     split_op(par, opr)
+
+    energy = calculate_energy(par, opr)
+    println("Energy is: ", energy)
 end
 
 main()