Developed Time Step Simulation

Here is a more developed form of the time step simulation. I've included the actual python script used to run it at the bottom of this post (if anybody knows a better way to post it let me know). Note that it is based on python 2, and requires the NumPy and Gnuplot python libraries to run. I've tried to document it well, but let me know if you have any questions.

Below are a few examples of what it can do. The first is a simulation of a perfect set of 1000 wells. The second is a simulation with the exact same parameters, but now with random gaussian disorder* of I/2 throughout, where I = 1 eV. The third has random gaussian disorder of I/3, and (as was requested by Mark) a linear adjustment to the wells, representing something like an electric field. Note that this last one also includes a new feature of a title with the disorder level and a countdown timer.

*disorder means the standard deviation of the gaussian distribution

Perfect Wells

Disorder I/2

Disorder I/3 and linear dependence

 

#simulation.py 

import numpy
from numpy import matrix
import math
import cmath
import Gnuplot
import random


#ESTABLISH CONSTANTS
#simulation constants (you probably don't need to change these)
I = 1.0 #eV
E0 = 16 #eV  ... this should not matter?
hbar = 6.582119 * ( 10**(-16) ) #eV*seconds
well_width = 2 #angstroms?
delta_t = hbar/(I) #seconds

#simulation parameters (change these at will)
disorder = I/3.0 #disorder of wells
time_steps = 2000 #steps of the simulation
m_size = 1000 #size of the matrix and vector
X0 = 500 #initial well that the wavefunction is localized in

#gnuplot parameters
output = 'output/simTest00.gif'
picture_count = 1000
delay = 2
#boundries of plot
xmin = 450
xmax = 550
ymin = 0
ymax = 0.2

eof = 'e' #end of file for gnuplot
#initialize gnuplot
GP = Gnuplot.Gnuplot()
GP( 'set term gif animate delay '+ str(delay) )
GP( 'set output \"' + output + '\" ' )
xr = '[' + str(xmin) + ':' + str(xmax) + ']' #x range
yr = '[' + str(ymin) + ':' + str(ymax) + ']' #y range
#END DEFINE CONSTANTS



#establish matrices
H     = numpy.asmatrix( numpy.zeros( shape=(m_size,m_size) ) ) #Hamiltonian
Ident = numpy.asmatrix( numpy.zeros( shape=(m_size,m_size) ) ) #Identity matrix
psi   = numpy.asmatrix( numpy.zeros( shape=(m_size,1) ) )      #wavefunction vector


#BEGIN DISORDER FUNCTIONS
#define any number of different disorder formats and choose them below

#produces random disorder with a gaussian distribution with
#a mean of 0 and a standard deviation of "disorder"
def random_gauss_disorder(disorder):
    return random.gauss(0,disorder)
#end random_gauss_disorder

#produces changes in the well in a linear progression
def linear_progression(well):
    multiplier = 1
    change = multiplier * well
    return change
#end linear_progession

#causes a sudden change in the wells at the edges
def cut_off(well):
    left_edge = 300
    right_edge = 800
    change = 0
    if (left_edge < well) and (well < right_edge):
        return 0
    else:
        return change
#end cut_off

#END DISORDER FUNCTIONS


#BEGIN MATRIX AND VECTOR INITIALIZATION
psi[X0] = 1
for i in range(0,m_size):
    #fill diagonal elements of H
    #THIS IS WHERE YOU PICK THE DISORDER FUNCTIONS (you can use several)  
    E = E0 + random_gauss_disorder(disorder) + linear_progression(i)
    H[i,i] = E
    #fill off-diagonal elements of H
    if not i == 0: H[i-1,i] = I
    if not i == m_size-1: H[i+1,i] = I
    #fill diagonal elements of Identity matrix
    Ident[i,i] = 1

M = Ident - ( 1j*(delta_t/hbar)*H ) # multiply by 1j for sqrt(-1)



print H
print
print M

print psi
print

#PERFORM TIME STEP ITERATIONS
picture_step = time_steps / picture_count
pic_counter  = 0

for step in range (0,time_steps):
    #time step iteration multiplication
    psi = M*psi

    #normalize vector
    con_tran = psi.transpose().conjugate()
    inner_product = (con_tran*psi)[0,0].real
    norm = 1 / math.sqrt(inner_product)
    psi = norm*psi

 
    if ( step == picture_step*pic_counter ):
        #find absolute squared amplitude of wave vector
        amp_list  = []
        for p in psi.tolist():
            sqr_real = math.pow(p[0].real,2)
            sqr_imag = math.pow(p[0].imag,2)
            abs_sqr = sqr_real + sqr_imag
            amp_list.append(abs_sqr)
     
     
        #turn amp_list into point format for gnuplot
        point_string_list = []
        for i in range( 0, len(amp_list) ):
            x = str(i)
            y = str( amp_list[i] )
         
            point = x + ' ' + y
            point_string_list.append(point)
     
        #plot data
        title = 'time = ' + str(time_steps-step) + ' ; disorder ' + str(disorder)
        GP('plot '+xr+yr+'\"-\" with boxes title \"' + title + '\"' )
        for point in point_string_list: GP(point)
        GP(eof)
     
        pic_counter+=1
        print step        #END TIME STEP ITERATIONS