Simulate Quantum Particles with Python and Qiskit
Code walkthrough on a classical and quantum computer.
View the complete code on my GitHub
Since their birth, quantum computers have been seen as a promising tool to be able to simulate what is happening at the quantum scale. In the status quo, we rely on computationally inefficient algorithms attempting to simulate quantum phenomena. Instead, why don't we just simulate quantum phenomena with the power of … quantum phenomena?
Quantum computers utilize particles on the quantum scale to store and process information. They have the unique ability to take on a continuum of states and generate interference between qubits. Utilizing these computers, simulating quantum mechanics becomes exponentially more efficient.
(For an in-depth refresher on the basics of quantum computers, refer to this article.)
Those who are familiar with basic quantum mechanics know that extremely small particles don’t behave in a predictable way. If I kick a soccer ball, I know exactly where it is going and where it is going to end up. At a quantum scale, the motion of small particles (i.e. electrons) is governed by, what is mathematically best described as, waves of probability. Formally, we call these “waves of probability” wave functions.
The graph below is a visual representation of a wave function. The amplitude of the wave relates to the probability of the particle being at that point in space.
This article provides the math and code behind using the split operator method to simulate the evolution of a particle’s wave function on a quantum computer. Feel free to skip around in the article to whatever section is most interesting to you:
1) Exploring the math behind the split operator method [click here]
2) Coding the algorithm on a classical computer (with python) [click here]
3) Coding the algirthm on a quantum computer (with python/qiskit) [click here]
The Math Behind the Split Operator Method
Erwin Shrödinger gave us the Shrödinger equation to describe an electron’s wavefunction based on external conditions.
In this differential equation, ℏ is the reduced Planck’s constant, Ĥ is a time-independent Hamiltonian, and ψ is the wavefunction of the particle (what we are solving for). A Hamiltonian represents the “total energy” of a system, including both its kinetic and potential energy.
We can integrate the equation on both sides to get the following:
We can now use this equation to find a particle’s wavefunction at time t based on its initial wavefunction and the hamiltonian of the system.
The term exp(-iĤt/ℏ) is used to calculate how the wave “propagates” over time, and is appropriately named the propagator.
As previously mentioned, Ĥ can be broken down into kinetic (T) and potential (V) energy:
Our high school physics knowledge covers the basics of breaking down kinetic energy in terms of momentum (p):
In an ideal world we can break down our propagator into kinetic and potential energy in one easy step:
Unfortunately, the two operators don’t commute. Concretely, exp(T+V) is not equal to exp(T)*exp(V). However, the Lie-Trotter product formula gives us an alternate way to split the two:
This essentially splits the propagator into an infinite amount of smaller steps alternating between the kinetic and potential propagators. The smaller the step, the more accurate our result will be. To reduce error, we can further split this by applying a half step of the potential propagator, a full step of the kinetic propagator, and the second half step of the potential propagator:
Substituting this into our integrated Shrödinger equation gives us:
Projecting this onto a coordinate basis ∣x⟩ gives us:
which can be simplified into our algorithm:
Note that we apply each operator on the RHS from right to left. The two integrals are just Fourier and inverse Fourier transforms to convert between the coordinate and momentum basis.
To translate the formula above, the algorithm will consist of 5 steps per time step/iteration:
- Apply a half step of the potential propagator to the initial state ψ
- Fourier transform into the momentum basis
- Apply a full step of the kinetic propagator on the momentum basis
- Inverse Fourier transform back into the coordinate basis
- Apply another half step of the potential propagator
*Note: I learned the math behind the algorithm from this video
Using Python to Code the Algorithm
Coding this on a classical computer is relatively straightforward. The equation above gives us our entire algorithm.
Here are all the libraries we are going to use. The first two lines are packages to better visualize+animate the wavefunction. Numpy allows us to work with arrays and makes the math a lot simpler to code.
from matplotlib import animation, rc, pyplot as pltfrom IPython.display import HTMLimport numpy as np
Next, let’s introduce some constants that we are going to be using in our code.
'''constants'''N=500 #number of grid pointshbar=1 #reduced plancks constantdt=0.5 #delta t, or the time stepm=1 #mass of particle in au
We then initialize spatial (x) and momentum (p) grids to apply the potential and kinetic propagators respectively. They are defined as follows:
np.linspace returns a NumPy array with N evenly spaced numbers between xMin and xMax, whereas np.arrange returns a NumPy array of all integers between -N/2 and N/2 with a default step size of 1.
'''spatial grid'''xMin=0xMax=500x=np.linspace(xMin, xMax, N)'''momentum grid'''p=np.arange(-N/2, N/2)*((2*np.pi*hbar)/(N*((xMax-xMin)/(N-1))))
Now that we have our constants and spatial/momentum grids, we can start our algorithm by creating the initial wavefunction of the particle.
The initial wavefunction is defined by the following equation:
where sigma can be interpreted as controlling the width of the initial wave, x0 as the initial position, and k0 as the initial momentum. Coding this equation is as simple as typing exactly what you see above.
'''wavefunction: guassian wavepacket'''k0=m/2 #initial momentumx0=xMax/2 #initial positionsigma=(xMax-xMin)/30 #widthpsi=(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-.5*((x-x0)/sigma)**2)*np.exp(1j*k0*x)
Graphing psi with matplotlib gives us the following:
Next, we have to define the propagators used in each iteration of the algorithm, beginning with the potential propagator. The potential energy of the system can be defined any way you like. To simulate a free particle set the potential energy to 0, and to create barriers set the potential energy to a certain value within the interval in space you would like the barrier to exist. The half step of the potential propagator is defined by the exponential derived in the first section of this article: exp(-iVt/2h).
'''potential operator+propogator'''Vhat=np.zeros_like(x) #free particleVhat[400:499]=1 #add well between x=[400,499] (optional)Vprop=np.exp(-1j*dt*Vhat/(2*hbar)) #half step
Graphing psi and V gives us the following:
Next, we define the kinetic propagator, given by the exponential derived in the first section: exp(-iTt/h) where T=p²/m.
'''kinetic operator+propogator'''That=p**2/(2*m)Tprop=np.exp(-1j*dt*That/hbar)
Finally, we define one time-step of our algorithm using the equation derived in the first section. A 1/2 step of the potential propagator, Fourier transform, a full step of the kinetic propagator, inverse Fourier transform, and finally, 1/2 step of the potential propagator is applied to the wavefunction psi.
np.fft.fft and np.fft.ifft apply a fast Fourier transform and inverse fast Fourier transform respectively.
'''function for one step of split operator method'''def splitOperator(psi,Vprop=Vprop,Tprop=Tprop): psi*=Vprop #1/2 step of potential propagator psi_p=np.fft.fft(psi) #shift to momentum basis psi_p*=Tprop #kinetic propagator psi=np.fft.ifft(psi_p) #shift back to spatial basis psi*=Vprop #1/2 step of potential propagator return psi
The algorithm is finished! Taking the initial wavefunction and updating it with the splitOperator() function for the desired number of iterations will simulate the evolution of the quantum state ψ. The following code animates the evolution of ψ in Python:
'''animate'''# First set up the figure, the axis, and the plot element we want to animatefig, ax = plt.subplots()plt.close()### parameters for plotax.set_xlim(( 0, 500))ax.set_ylim((-0.001, 0.001))line, = ax.plot([], [], lw=2)lineV, = ax.plot([], [], lw=2)# initialization function: plot the background of each framedef init(): line.set_data([], []) lineV.set_data([], []) return (line, lineV,)def animate(i): global psi psi=splitOperator(psi) line.set_data(x, np.conj(psi)*psi) lineV.set_data(x, Vhat) return(line, lineV,)anim = animation.FuncAnimation(fig, animate, init_func=init,frames=400, interval=10, blit=True)# Note: below is the part which makes it work on Colabrc('animation', html='jshtml')anim
Using Qiskit to Code the Algorithm on a Quantum Computer
The algorithm on a quantum computer is similar to the algorithm we just coded, except we do not split the potential propagator into two half steps. The steps are listed below:
- Initialize the wavefunction discretized into 2^N gridpoints, where N is the number of qubits (excluding ancillary) used in the circuit.
- Apply the potential propagator exp(-iVt).
- Apply a Quantum Fourier Transform (QFT)
- Apply the kinetic propagator exp(-iTt)
- Apply an inverse QFT
- Repeat steps 2–5 for the desired number of iterations
- Measure the circuit
Hardware on current quantum computers is extremely limited. Thus, we can’t casually create 500 grid points (using 9 qubits) in our spatial grid without incurring a large amount of error. For this code, we are using 7 qubits, allowing for 2⁶ (64) grid points and 1 ancillary qubit.
Here are all the libraries we import to run the code. Qiskit is a library that allows us to code on real or simulated quantum computers. Usually, pylatexenc isn’t required as an import, however, on Google Colab I needed to write a separate import statement to be able to visualize my quantum circuits. (specifically running circuit.draw(‘mpl’) where the circuit is just any quantum circuit object)
from qiskit import *import numpy as npfrom matplotlib import pyplot as pltfrom qiskit.visualization import plot_histogramfrom qiskit.circuit.library import QFTimport pylatexenc
Next, we define some constants:
'''constants'''n=6 #number of qubitsN=2**n #grid pointsphi=0.1 #time stepiterations=50
and the potential energy:
'''potential'''#V=0 free particleV=1.75
Note controlling the location of any wells/barriers is done by controlling what gates you apply the potential propagator to.
To apply exp(i*V*phi) to a qubit, we apply a phase gate parameterized by the value V*phi.
For example, to create a well in the second half of the spatial grid apply the potential propagator on the last qubit of the circuit. This is essentially applied to all values between 100000–111111, or x[32:64].
*Note that the first qubit is the ancillary qubit
To create a barrier from 001110 to 001111, we can apply a series of Toffoli and a controlled phase gate to apply the potential propagator on the 3rd, 4th, and 5th qubits.
Next, we must determine how to represent the kinetic propagator, exp(-iTt), where T=p²/m. We can write p as follows:
Thus the whole propagator exp(-ip²t) can be written as follows:
where Z is a phase gate that applies the following operation onto qubit j:
By expanding the squared term in the previous equation it is shown that Z operators are acting on all individual qubits (excluding ancillary) and commuting Z operators acting on all qubit pairs as well. With a phi value of 0.1 and a 7 qubit circuit (q0 is an ancillary qubit), the kinetic propagator looks like this:
The next step is to code a single time step of the split operator algorithm using the propagators we defined above:
qubits=QuantumRegister(n+1)splitOp=QuantumCircuit(qubits,name="splitOp")
'''potential operator'''#adding a barrier from 001110 to 001111splitOp.ccx(3,4,0)splitOp.cp(V*phi,0,5)splitOp.ccx(3,4,0)splitOp.barrier()
'''QFT'''splitOp.append(QFT(n),qubits[1:])
'''momentum operator'''for i in range(n): splitOp.p(phi/2**(n-3+i), n-i)for i in range(n): for j in range(i+1,n): splitOp.cx(n-j,0) splitOp.cx(n-i,0) splitOp.p(phi*(2**(2-i-j)),0) splitOp.cx(n-i,0) splitOp.cx(n-j,0)
'''inverse QFT'''splitOp.append(QFT(n,inverse=True),qubits[1:]) #inverse QFTsplitOp.barrier()
#splitOp.draw('mpl')
We then build our final circuit based on the number of iterations defined at the beginning of the code. First, we initialize the gaussian wavepacket onto our qubits and then add the splitOp circuit for each time step. Each qubit (minus the ancillary) is measured onto a classical bit. The circuit is run on Qiskit’s qasm_simulator and stored in the answer dictionary:
q=QuantumRegister(n+1)c=ClassicalRegister(n)circuit=QuantumCircuit(q,c)
'''initialization'''circuit.initialize(psi,q[1:])circuit.barrier()
'''time steps'''for i in range(iterations): circuit.append(splitOp,q[:])
'''measurement'''for i in range(1,n+1): circuit.measure(i,i-1)
#circuit.draw('mpl')answer = execute(circuit, backend=BasicAer.get_backend('qasm_simulator'), shots=5000).result().get_counts()
The following code plots the results on a histogram:
for a in range(N): s = bin(a)[2:] while len(s)!=6: s='0'+s if s not in answer.keys(): answer[s] = 0
plot_histogram(answer, figsize = [20,10])
Above we can see a particle tunneling through a barrier from 001110 to 001111.
What’s Next?
After learning the basics of how to simulate one particle moving with the split operator method, the next step is looking at interactions between particles and chemical dynamics. Research in particle simulation has applications in drug development, condensed matter physics, material science, etc.
Thank you for reading my article! If you enjoyed reading, press the clap button below. Here are links to connect with me and see what I am up to:
Github
Email: lohitpotnuru@gmail.com
References
https://www.youtube.com/watch?v=o96K8fkOrG8&t=327s
https://www.youtube.com/watch?v=7Bsx_4hOU3g&t=4470s
https://fount.aucegypt.edu/cgi/viewcontent.cgi?article=1020&context=capstone
More content at plainenglish.io. Sign up for our free weekly newsletter. Get exclusive access to writing opportunities and advice in our community Discord.
