Initial Conditions

Our simulation evolves the Schrödinger equation in time, but without initial conditions it has nothing to evolve, and is a very boring simulation. Let's spice things up a bit and fire in a free particle from the left. Later, we will put the variations in the potential, $V\left(x\right)$ , to the right so that the particle will encounter and interact with the potential.

The free particle Schrödinger equation $$i\hslash \frac{\partial }{\partial t}\Psi \left(x,t\right)=\frac{-{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}\Psi \left(x,t\right)$$ has solutions representing a localized free particle as a Gaussian wave packet. A Gaussian wave packet centered around $x=x_0$ with momentum $p_0$ is¹ $$\Psi \left(x,0\right)=\frac{1}{{\left(2\pi {\mathrm{\sigma }}^2\right)}^{\frac{1}{4}}}{e}^{-\frac{{(x-x_0)}^2}{4{\mathrm{\sigma }}^2}}{e}^{i{p}_0\mathrm{x}}$$

The most interesting part is the momentum exponential, where $p=\hslash k$ , but in our code $\hslash =1$ , so we make the expansion $$e^{i{p}_0\mathrm{x}}=\cos kx+i\sin kx$$ This is the only term with an imaginary part, so it determines the division between real and imaginary parts of the wave function. We represent the wave function as an array<vec2f>, with the 0 element being the real part, and the 1 element being the imaginary part. We set components of our wave function to this phase value, then multiply the vec2 by the magnitude of the wave function. The shader language, WGSL, takes care of multiplying each of the vec2 components by this magnitude factor.

There are no derivatives in this expression, so we need only evaluate if for each position on our grid. Luckily that is exactly what we have been doing.

      const PI : f32 = 3.1415926535897932384626433832795;

      struct Parameters {
        dt: f32,          // The time step size
        xResolution: u32, // The number of points along the x-axis, the number of elements in the array.
        length: f32       // The full length for our simulation
      }

      struct WavePacketParameters {
        x0: f32, // The center of the wave packet.
        w: f32,  // The width of the wave packet, in GeV^-1.
        k: f32   // The wave number, from particle momentum.
      }

      // group 0, parameters that never change over the course of a simulation
      @group(0) @binding(0) var<storage, read> parameters: Parameters;

      // Group 1, parameters describing the wave packet
      @group(1) @binding(0) var<storage, read> wavePacketParameters : WavePacketParameters;

      // Initial wave function at t=0.
      @group(2) @binding(0) var<storage, read_write> waveFunction : array<vec2f>;

      fn computePsi(globalID: u32) -> vec2f
      {
        let x = f32(globalID) * parameters.length / f32(parameters.xResolution);
        let alpha = wavePacketParameters.w * wavePacketParameters.w;
        let deltaX = x - wavePacketParameters.x0;
        // Normalization constant http://quantummechanics.ucsd.edu/ph130a/130_notes/node80.html
        let a = pow(2.0/(PI*alpha), 0.25);
        let gaussian = a*exp(-deltaX*deltaX/alpha);
        let phase = vec2(cos(wavePacketParameters.k*x), sin(wavePacketParameters.k*x));

        return gaussian*phase;
      }

      @compute @workgroup_size(64)
      fn computeInitialValues(@builtin(global_invocation_id) global_id : vec3u)
      {
        let index = global_id.x;
        // Skip invocations when work groups exceed the actual problem size
        if (index >= parameters.xResolution) {
          return;
        }
        waveFunction[index] = computePsi(index);
      }
      

There are some interesting things we can observe in even this initial simulation. The wave packet is not a sharp peak, but is spread out in space. Further, the wave packet spreads out as it propagates, so it is wider and shorter at the end of the simulation. That is uncertainty in the location of a particle naturally increases over time even with no interactions.

An important technical aspect is the performance of this simulation. It seems extremely slow, especially when compared with an earlier WebGL implementation. The next section looks at some performance measurements and tuning.

• Back
• Next