Quick Start Guide¶
Simulation Set up¶
Start by importing the module:
import trottersuzuki as ts
To set up the simulation of a quantum system, we need only a few lines of code.
First of all, we create the lattice over which the physical system is
defined. All information about the discretized space is collected in a
single object. Say we want a squared lattice of 300x300 nodes, with a
physical area of 20x20, then we have to specify these in the constructor of the Lattice
class:
grid = ts.Lattice2D(300, 20.)
The object grid
defines the geometry of the system and it
will be used throughout the simulations. Note that the origin of the lattice is at its centre and the lattice points are scaled to the physical locations.
The physics of the problem is described by the Hamiltonian. A single object is going to store all the information regarding the Hamiltonian. The module is able to deal with two physical models: Gross-Pitaevskii equation of a single or two-component wave function, namely (in units \(\hbar=1\)):
being, for Cartesian coordinates
while for cylindrical coordinates
and \(\psi(t) = \psi_t(x,y)\) for the single component wave function, or
where
and \(\psi(t) = \begin{bmatrix} \psi_1(t) \\ \psi_2(t) \end{bmatrix}\), for the two component wave function.
First we define the object for the external potential \(V(x,y)\). A general external potential function can be defined by a Python function, for instance, the harmonic potential can be defined as follows:
def harmonic_potential(x,y):
return 0.5 * (x**2 + y**2)
Now we create the external potential object using the Potential
class and then we initialize it with the function above:
potential = ts.Potential(grid) # Create the potential object
potential.init_potential(harmonic_potential) # Initialize it using a python function
Note that the module provides a quick way to define the harmonic potential, as it is fequently used:
omegax = omegay = 1.
harmonicpotential = ts.HarmonicPotential(grid, omegax, omegay)
We are ready to create the Hamiltonian
object. For the sake of simplicity, let us create the Hamiltonian of the harmonic oscillator:
particle_mass = 1. # Mass of the particle
hamiltonian = ts.Hamiltonian(grid, potential, particle_mass) # Create the Hamiltonian object
The quantum state is created by the State
class; it resembles the way the potential is defined. Here we create the ground state of the
harmonic oscillator:
import numpy as np # Import the module numpy for the exponential and sqrt functions
def state_wave_function(x,y): # Wave function
return np.exp(-0.5*(x**2 + y**2)) / np.sqrt(np.pi)
state = ts.State(grid) # Create the quantum state
state.init_state(state_wave_function) # Initialize the state
The module provides several predefined quantum states as well. In this
case, we could have used the GaussianState
class:
omega = 1.
gaussianstate = ts.GaussianState(grid, omega) # Create a quantum state whose wave function is Gaussian-like
We are left with the creation of the last object: the Solver
class gathers all the objects we defined so far and it is used to perform the evolution and analyze the expectation values:
delta_t = 1e-3 # Physical time of a single iteration
solver = ts.Solver(grid, state, hamiltonian, delta_t) # Creating the solver object
If you want to solve the problem on the GPU, request the matching kernel:
delta_t = 1e-3 # Physical time of a single iteration
solver = ts.Solver(grid, state, hamiltonian, delta_t, kernel_type="gpu")
Note that not all functionality is available in the GPU kernel.
Finally we can perform both real-time and imaginary-time evolution using
the method evolve
:
iterations = 100 # Number of iterations to be performed
solver.evolve(iterations, True) # Perform imaginary-time evolution
solver.evolve(iterations) # Perform real-time evolution
Analysis¶
The classes we have seen so far implement several members useful to analyze the system (see the function reference section for a complete list).
Expectation values¶
The solver class provides members for the energy calculations. For
instance, the total energy can be calculated using the
get_total_energy
member. We expect it to be \(1\)
(\(\hbar =1\)), and indeed we get the right result up to a small
error which depends on the lattice approximation:
tot_energy = solver.get_total_energy()
print(tot_energy)
1.00146456951
The expected values of the \(X\), \(Y\), \(P_x\),
\(P_y\) operators are calculated using the members in the State
class
mean_x = state.get_mean_x() # Get the expected value of X operator
print(mean_x)
1.39431975344e-14
Norm of the state¶
The squared norm of the state can be calculated by means of both
State
and Solver
classes
snorm = state.get_squared_norm()
print(snorm)
1.0
Particle density and Phase¶
Very often one is interested in the phase and particle density of the
state. Two members of State
class provide these features
density = state.get_particle_density() # Return a numpy matrix of the particle density
phase = state.get_phase() # Return a numpy matrix of the phase
Imprinting¶
The member imprint
, in the State
class, applies the following transformation to the state:
being \(f(x,y)\) a general complex-valued function. This comes in handy when we want to imprint, for instance, vortices or solitons:
def vortex(x, y): # Function defining a vortex
z = x + 1j*y
angle = np.angle(z)
return np.exp(1j * angle)
state.imprint(vortex) # Imprint the vortex on the state
File Input and Output¶
write_to_files
and loadtxt
members, in State
class, provide
a simple way to handle file I/O. The former writes the wave function
arranged as a complex matrix, in a plain text; the latter loads the wave
function from a file to the state
object. The following code
provides an example:
state.write_to_file("file_name") # Write the wave function to a file
state2 = ts.State(grid) # Create a new state
state2.loadtxt("file_name") # Load the wave function from the file
For a complete list of methods see the function reference.