2+2 (50-5*6)/4 7/3 7/3. 7/float(3) sqrt(81) import math math.sqrt(81) width = 20 length = 30 area = length*width area volume depth = 10 volume = area*depth volume return = 0 'Hello, World!' "Hello, World!" "He's a Rebel" 'She asked, "How are you today?"' greeting = "Hello, World!" print greeting print "The area is ",area statement = "Hello," + "World!" print statement statement = "Hello, " + "World!" print statement print "This " + "is " + "a " + "longer " + "statement." days_of_the_week = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"] days_of_the_week[2] days_of_the_week[-1] languages = ["Fortran","C","C++"] languages.append("Python") print languages range(10) range(2,8) evens = range(0,20,2) evens evens[3] ["Today",7,99.3,""] help(len) len(evens) for day in days_of_the_week: print day for day in days_of_the_week: statement = "Today is " + day print statement for i in range(20): print "The square of ",i," is ",i*i for letter in "Sunday": print letter days_of_the_week[0] days_of_the_week[0:2] days_of_the_week[:2] days_of_the_week[-2:] workdays = days_of_the_week[1:6] print workdays day = "Sunday" abbreviation = day[:3] print abbreviation numbers = range(0,40) evens = numbers[2::2] evens if day == "Sunday": print "Sleep in" else: print "Go to work" day == "Sunday" 1 == 2 50 == 2*25 3 < 3.14159 1 == 1.0 1 != 0 1 <= 2 1 >= 1 1 is 1.0 [1,2,3] == [1,2,4] [1,2,3] < [1,2,4] hours = 5 0 < hours < 24 if day == "Sunday": print "Sleep in" elif day == "Saturday": print "Do chores" else: print "Go to work" for day in days_of_the_week: statement = "Today is " + day print statement if day == "Sunday": print " Sleep in" elif day == "Saturday": print " Do chores" else: print " Go to work" bool(1) bool(0) bool(["This "," is "," a "," list"]) n = 10 sequence = [0,1] for i in range(2,n): # This is going to be a problem if we ever set n <= 2! sequence.append(sequence[i-1]+sequence[i-2]) print sequence def fibonacci(sequence_length): "Return the Fibonacci sequence of length *sequence_length*" sequence = [0,1] if sequence_length < 1: print "Fibonacci sequence only defined for length 1 or greater" return if 0 < sequence_length < 3: return sequence[:sequence_length] for i in range(2,sequence_length): sequence.append(sequence[i-1]+sequence[i-2]) return sequence fibonacci(2) fibonacci(12) help(fibonacci) from math import factorial help(factorial) factorial(20) def fact(n): if n <= 0: return 1 return n*fact(n-1) fact(20) t = (1,2,'hi',9.0) t t[1] t.append(7) t[1]=77 ('Bob',0.0,21.0) positions = [ ('Bob',0.0,21.0), ('Cat',2.5,13.1), ('Dog',33.0,1.2) ] def minmax(objects): minx = 1e20 # These are set to really big numbers miny = 1e20 for obj in objects: name,x,y = obj if x < minx: minx = x if y < miny: miny = y return minx,miny x,y = minmax(positions) print x,y x,y = 1,2 y,x = x,y x,y mylist = [1,2,9,21] ages = {"Rick": 46, "Bob": 86, "Fred": 21} print "Rick's age is ",ages["Rick"] dict(Rick=46,Bob=86,Fred=20) len(t) len(ages) fibs = fibonacci(10) facts = [] for i in range(10): facts.append(factorial(i)) figsize(8,6) plot(facts,label="factorial") plot(fibs,label="Fibonacci") xlabel("n") legend() semilogy(facts,label="factorial") semilogy(fibs,label="Fibonacci") xlabel("n") legend() import this array([1,2,3,4,5,6]) array([1,2,3,4,5,6],'d') array([1,2,3,4,5,6],'D') array([1,2,3,4,5,6],'i') array([[0,1],[1,0]],'d') zeros((3,3),'d') zeros(3,'d') zeros((1,3),'d') zeros((3,1),'d') identity(4,'d') linspace(0,1) linspace(0,1,11) x = linspace(0,2*pi) sin(x) plot(x,sin(x)) 0.125*identity(3,'d') identity(2,'d') + array([[1,1],[1,2]]) identity(2)*ones((2,2)) dot(identity(2),ones((2,2))) v = array([3,4],'d') sqrt(dot(v,v)) m = array([[1,2],[3,4]]) m.T diag([1,2,3,4,5]) A = array([[1,1,1],[0,2,5],[2,5,-1]]) b = array([6,-4,27]) solve(A,b) A = array([[13,-4],[-4,7]],'d') eigvalsh(A) eigh(A) def nderiv(y,x): "Finite difference derivative of the function f" n = len(y) d = zeros(n,'d') # assume double # Use centered differences for the interior points, one-sided differences for the ends for i in range(1,n-1): d[i] = (y[i+1]-y[i-1])/(x[i+1]-x[i-1]) d[0] = (y[1]-y[0])/(x[1]-x[0]) d[n-1] = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]) return d x = linspace(0,2*pi) dsin = nderiv(sin(x),x) plot(x,dsin,label='numerical') plot(x,cos(x),label='analytical') title("Comparison of numerical and analytical derivatives of sin(x)") legend() def Laplacian(x): h = x[1]-x[0] # assume uniformly spaced points n = len(x) M = -2*identity(n,'d') for i in range(1,n): M[i,i-1] = M[i-1,i] = 1 return M/h**2 x = linspace(-3,3) m = 1.0 ohm = 1.0 T = (-0.5/m)*Laplacian(x) V = 0.5*(ohm**2)*(x**2) H = T + diag(V) E,U = eigh(H) h = x[1]-x[0] # Plot the Harmonic potential plot(x,V,color='k') for i in range(4): # For each of the first few solutions, plot the energy level: axhline(y=E[i],color='k',ls=":") # as well as the eigenfunction, displaced by the energy level so they don't # all pile up on each other: plot(x,-U[:,i]/sqrt(h)+E[i]) title("Eigenfunctions of the Quantum Harmonic Oscillator") xlabel("Displacement (bohr)") ylabel("Energy (hartree)") from numpy.polynomial.hermite import Hermite def ho_evec(x,n,m,ohm): vec = [0]*9 vec[n] = 1 Hn = Hermite(vec) return (1/sqrt(2**n*factorial(n)))*pow(m*ohm/pi,0.25)*exp(-0.5*m*ohm*x**2)*Hn(x*sqrt(m*ohm)) plot(x,ho_evec(x,0,1,1),label="Analytic") plot(x,-U[:,0]/sqrt(h),label="Numeric") xlabel('x (bohr)') ylabel(r'$\psi(x)$') title("Comparison of numeric and analytic solutions to the Harmonic Oscillator") legend() phase_correction = [-1,1,1,-1,-1,1] for i in range(6): subplot(2,3,i+1) plot(x,ho_evec(x,i,1,1),label="Analytic") plot(x,phase_correction[i]*U[:,i]/sqrt(h),label="Numeric") from scipy.special import airy,jn,eval_chebyt,eval_legendre subplot(2,2,1) x = linspace(-1,1) Ai,Aip,Bi,Bip = airy(x) plot(x,Ai) plot(x,Aip) plot(x,Bi) plot(x,Bip) title("Airy functions") subplot(2,2,2) x = linspace(0,10) for i in range(4): plot(x,jn(i,x)) title("Bessel functions") subplot(2,2,3) x = linspace(-1,1) for i in range(6): plot(x,eval_chebyt(i,x)) title("Chebyshev polynomials of the first kind") subplot(2,2,4) x = linspace(-1,1) for i in range(6): plot(x,eval_legendre(i,x)) title("Legendre polynomials") raw_data = """\ 3.1905781584582433,0.028208609537968457 4.346895074946466,0.007160804747670053 5.374732334047101,0.0046962988461934805 8.201284796573875,0.0004614473299618756 10.899357601713055,0.00005038370219939726 16.295503211991434,4.377451812785309e-7 21.82012847965739,3.0799922117601088e-9 32.48394004282656,1.524776208284536e-13 43.53319057815846,5.5012073588707224e-18""" data = [] for line in raw_data.splitlines(): words = line.split(',') data.append(map(float,words)) data = array(data) title("Raw Data") xlabel("Distance") plot(data[:,0],data[:,1],'bo') title("Raw Data") xlabel("Distance") semilogy(data[:,0],data[:,1],'bo') params = polyfit(data[:,0],log(data[:,1]),1) a = params[0] A = exp(params[1]) x = linspace(1,45) title("Raw Data") xlabel("Distance") semilogy(data[:,0],data[:,1],'bo') semilogy(x,A*exp(a*x),'b-') gauss_data = """\ -0.9902286902286903,1.4065274110372852e-19 -0.7566104566104566,2.2504438576596563e-18 -0.5117810117810118,1.9459459459459454 -0.31887271887271884,10.621621621621626 -0.250997150997151,15.891891891891893 -0.1463309463309464,23.756756756756754 -0.07267267267267263,28.135135135135133 -0.04426734426734419,29.02702702702703 -0.0015939015939017698,29.675675675675677 0.04689304689304685,29.10810810810811 0.0840994840994842,27.324324324324326 0.1700546700546699,22.216216216216214 0.370878570878571,7.540540540540545 0.5338338338338338,1.621621621621618 0.722014322014322,0.08108108108108068 0.9926849926849926,-0.08108108108108646""" data = [] for line in gauss_data.splitlines(): words = line.split(',') data.append(map(float,words)) data = array(data) plot(data[:,0],data[:,1],'bo') def gauss(x,A,a): return A*exp(a*x**2) from scipy.optimize import curve_fit params,conv = curve_fit(gauss,data[:,0],data[:,1]) x = linspace(-1,1) plot(data[:,0],data[:,1],'bo') A,a = params plot(x,gauss(x,A,a),'b-') from random import random rands = [] for i in range(100): rands.append(random()) plot(rands) from random import gauss grands = [] for i in range(100): grands.append(gauss(0,1)) plot(grands) plot(rand(100)) npts = 5000 xs = 2*rand(npts)-1 ys = 2*rand(npts)-1 r = xs**2+ys**2 ninside = (r<1).sum() figsize(6,6) # make the figure square title("Approximation to pi = %f" % (4*ninside/float(npts))) plot(xs[r<1],ys[r<1],'b.') plot(xs[r>1],ys[r>1],'r.') figsize(8,6) # change the figsize back to 4x3 for the rest of the notebook n = 100 total = 0 for k in range(n): total += pow(-1,k)/(2*k+1.0) print 4*total from numpy import sqrt def f(x): return exp(-x) x = linspace(0,10) plot(x,exp(-x)) from scipy.integrate import quad quad(f,0,inf) from scipy.fftpack import fft,fftfreq npts = 4000 nplot = npts/10 t = linspace(0,120,npts) def acc(t): return 10*sin(2*pi*2.0*t) + 5*sin(2*pi*8.0*t) + 2*rand(npts) signal = acc(t) FFT = abs(fft(signal)) freqs = fftfreq(npts, t[1]-t[0]) subplot(211) plot(t[:nplot], signal[:nplot]) subplot(212) plot(freqs,20*log10(FFT),',') show() myoutput = """\ @ Step Energy Delta E Gmax Grms Xrms Xmax Walltime @ ---- ---------------- -------- -------- -------- -------- -------- -------- @ 0 -6095.12544083 0.0D+00 0.03686 0.00936 0.00000 0.00000 1391.5 @ 1 -6095.25762870 -1.3D-01 0.00732 0.00168 0.32456 0.84140 10468.0 @ 2 -6095.26325979 -5.6D-03 0.00233 0.00056 0.06294 0.14009 11963.5 @ 3 -6095.26428124 -1.0D-03 0.00109 0.00024 0.03245 0.10269 13331.9 @ 4 -6095.26463203 -3.5D-04 0.00057 0.00013 0.02737 0.09112 14710.8 @ 5 -6095.26477615 -1.4D-04 0.00043 0.00009 0.02259 0.08615 20211.1 @ 6 -6095.26482624 -5.0D-05 0.00015 0.00002 0.00831 0.03147 21726.1 @ 7 -6095.26483584 -9.6D-06 0.00021 0.00004 0.01473 0.05265 24890.5 @ 8 -6095.26484405 -8.2D-06 0.00005 0.00001 0.00555 0.01929 26448.7 @ 9 -6095.26484599 -1.9D-06 0.00003 0.00001 0.00164 0.00564 27258.1 @ 10 -6095.26484676 -7.7D-07 0.00003 0.00001 0.00161 0.00553 28155.3 @ 11 -6095.26484693 -1.8D-07 0.00002 0.00000 0.00054 0.00151 28981.7 @ 11 -6095.26484693 -1.8D-07 0.00002 0.00000 0.00054 0.00151 28981.7""" lines = myoutput.splitlines() lines for line in lines[2:]: # do something with each line words = line.split() import string help(string.split) lines[2].split() for line in lines[2:]: # do something with each line words = line.split() energy = words[2] gmax = words[4] time = words[8] print energy,gmax,time data = [] for line in lines[2:]: # do something with each line words = line.split() energy = float(words[2]) gmax = float(words[4]) time = float(words[8]) data.append((energy,gmax,time)) data = array(data) plot(data[:,0]) xlabel('step') ylabel('Energy (hartrees)') title('Convergence of NWChem geometry optimization for Si cluster') csv = """\ -6095.12544083, 0.03686, 1391.5 -6095.25762870, 0.00732, 10468.0 -6095.26325979, 0.00233, 11963.5 -6095.26428124, 0.00109, 13331.9 -6095.26463203, 0.00057, 14710.8 -6095.26477615, 0.00043, 20211.1 -6095.26482624, 0.00015, 21726.1 -6095.26483584, 0.00021, 24890.5 -6095.26484405, 0.00005, 26448.7 -6095.26484599, 0.00003, 27258.1 -6095.26484676, 0.00003, 28155.3 -6095.26484693, 0.00002, 28981.7 -6095.26484693, 0.00002, 28981.7""" data = [] for line in csv.splitlines(): words = line.split(',') data.append(map(float,words)) data = array(data) help(map) plot(data[:,0]) xlabel('step') ylabel('Energy (hartrees)') title('Convergence of NWChem geometry optimization for Si cluster') energies = data[:,0] minE = min(energies) energies_eV = 27.211*(energies-minE) plot(energies_eV) xlabel('step') ylabel('Energy (eV)') title('Convergence of NWChem geometry optimization for Si cluster') lines = """\ ---------------------------------------- | WALL | 0.45 | 443.61 | ---------------------------------------- @ Step Energy Delta E Gmax Grms Xrms Xmax Walltime @ ---- ---------------- -------- -------- -------- -------- -------- -------- @ 0 -6095.12544083 0.0D+00 0.03686 0.00936 0.00000 0.00000 1391.5 ok ok Z-matrix (autoz) -------- """.splitlines() for line in lines: if line.startswith('@'): print line print "I have 3 errands to run" "I have 3 errands to run" a,b,c = 1,2,3 print "The variables are ",1,2,3 print "Pi as a decimal = %d" % pi print "Pi as a float = %f" % pi print "Pi with 4 decimal places = %.4f" % pi print "Pi with overall fixed length of 10 spaces, with 6 decimal places = %10.6f" % pi print "Pi as in exponential format = %e" % pi print "The variables specified earlier are %d, %d, and %d" % (a,b,c) form_letter = """\ %s Dear %s, We regret to inform you that your product did not ship today due to %s. We hope to remedy this as soon as possible. From, Your Supplier """ print form_letter % ("July 1, 2013","Valued Customer Bob","alien attack") form_letter = """\ %(date)s Dear %(customer)s, We regret to inform you that your product did not ship today due to %(lame_excuse)s. We hope to remedy this as soon as possible. From, Your Supplier """ print form_letter % {"date" : "July 1, 2013","customer":"Valued Customer Bob","lame_excuse":"alien attack"} nwchem_format = """ start %(jobname)s title "%(thetitle)s" charge %(charge)d geometry units angstroms print xyz autosym %(geometry)s end basis * library 6-31G** end dft xc %(dft_functional)s mult %(multiplicity)d end task dft %(jobtype)s """ oxygen_xy_coords = [(0,0),(0,0.1),(0.1,0),(0.1,0.1)] charge = 0 multiplicity = 1 dft_functional = "b3lyp" jobtype = "optimize" geometry_template = """\ O %f %f 0.0 H 0.0 1.0 0.0 H 1.0 0.0 0.0""" for i,xy in enumerate(oxygen_xy_coords): thetitle = "Water run #%d" % i jobname = "h2o-%d" % i geometry = geometry_template % xy print "---------" print nwchem_format % dict(thetitle=thetitle,charge=charge,jobname=jobname,jobtype=jobtype, geometry=geometry,dft_functional=dft_functional,multiplicity=multiplicity) def my_enumerate(seq): l = [] for i in range(len(seq)): l.append((i,seq[i])) return l my_enumerate(oxygen_xy_coords) linspace(0,1) linspace(0,1,5) linspace(0,1,5,endpoint=False) def my_linspace(start,end): npoints = 50 v = [] d = (end-start)/float(npoints-1) for i in range(npoints): v.append(start + i*d) return v my_linspace(0,1) def my_linspace(start,end,npoints = 50): v = [] d = (end-start)/float(npoints-1) for i in range(npoints): v.append(start + i*d) return v my_linspace(0,1) my_linspace(0,1,5) def my_linspace(start,end,npoints=50,**kwargs): endpoint = kwargs.get('endpoint',True) v = [] if endpoint: d = (end-start)/float(npoints-1) else: d = (end-start)/float(npoints) for i in range(npoints): v.append(start + i*d) return v my_linspace(0,1,5,endpoint=False) def my_range(*args): start = 0 step = 1 if len(args) == 1: end = args[0] elif len(args) == 2: start,end = args elif len(args) == 3: start,end,step = args else: raise Exception("Unable to parse arguments") v = [] value = start while True: v.append(value) value += step if value > end: break return v my_range() evens1 = [2*i for i in range(10)] print evens1 odds = [i for i in range(20) if i%2==1] odds def evens_below(n): for i in xrange(n): if i%2 == 0: yield i return for i in evens_below(9): print i list(evens_below(9)) evens_gen = (i for i in xrange(9) if i%2==0) for i in evens_gen: print i def gauss(x,A,a,x0): return A*exp(-a*(x-x0)**2) def gauss_maker(A,a,x0): def f(x): return A*exp(-a*(x-x0)**2) return f x = linspace(0,1) g = gauss_maker(99.0,1.0,0.5) plot(x,g(x)) # Data in a json format: json_data = """\ { "a": [1,2,3], "b": [4,5,6], "greeting" : "Hello" }""" import json json.loads(json_data) json.dumps({"a":[1,2,3],"b":[9,10,11],"greeting":"Hola"}) from operator import add, mul add(1,2) mul(3,4) def doubler(x): return 2*x doubler(17) lambda x: 2*x another_doubler = lambda x: 2*x another_doubler(19) map(float,'1 2 3 4 5'.split()) sum([1,2,3,4,5]) def prod(l): return reduce(mul,l) prod([1,2,3,4,5]) mystring = "Hi there" mystring.split() mystring.startswith('Hi') len(mystring) class Schrod1d: """\ Schrod1d: Solver for the one-dimensional Schrodinger equation. """ def __init__(self,V,start=0,end=1,npts=50,**kwargs): m = kwargs.get('m',1.0) self.x = linspace(start,end,npts) self.Vx = V(self.x) self.H = (-0.5/m)*self.laplacian() + diag(self.Vx) return def plot(self,*args,**kwargs): titlestring = kwargs.get('titlestring',"Eigenfunctions of the 1d Potential") xstring = kwargs.get('xstring',"Displacement (bohr)") ystring = kwargs.get('ystring',"Energy (hartree)") if not args: args = [3] x = self.x E,U = eigh(self.H) h = x[1]-x[0] # Plot the Potential plot(x,self.Vx,color='k') for i in range(*args): # For each of the first few solutions, plot the energy level: axhline(y=E[i],color='k',ls=":") # as well as the eigenfunction, displaced by the energy level so they don't # all pile up on each other: plot(x,U[:,i]/sqrt(h)+E[i]) title(titlestring) xlabel(xstring) ylabel(ystring) return def laplacian(self): x = self.x h = x[1]-x[0] # assume uniformly spaced points n = len(x) M = -2*identity(n,'d') for i in range(1,n): M[i,i-1] = M[i-1,i] = 1 return M/h**2 square_well = Schrod1d(lambda x: 0*x,m=10) square_well.plot(4,titlestring="Square Well Potential") ho = Schrod1d(lambda x: x**2,start=-3,end=3) ho.plot(6,titlestring="Harmonic Oscillator") def finite_well(x,V_left=1,V_well=0,V_right=1,d_left=10,d_well=10,d_right=10): V = zeros(x.size,'d') for i in range(x.size): if x[i] < d_left: V[i] = V_left elif x[i] > (d_left+d_well): V[i] = V_right else: V[i] = V_well return V fw = Schrod1d(finite_well,start=0,end=30,npts=100) fw.plot() def triangular(x,F=30): return F*x tw = Schrod1d(triangular,m=10) tw.plot() def tri_finite(x): return finite_well(x)+triangular(x,F=0.025) tfw = Schrod1d(tri_finite,start=0,end=30,npts=100) tfw.plot() %timeit factorial(20) %timeit fact(20) def evens(n): "Return a list of even numbers below n" l = [] for x in range(n): if x % 2 == 0: l.append(x) return l import cProfile cProfile.run('evens(100000)') def evens2(n): "Return a list of even numbers below n" return [x for x in range(n) if x % 2 == 0] import cProfile cProfile.run('evens2(100000)') def evens3(n): "Return a list of even numbers below n" return [x for x in xrange(n) if x % 2 == 0] import cProfile cProfile.run('evens3(100000)') def primes(n): """\ From python cookbook, returns a list of prime numbers from 2 to < n >>> primes(2) [2] >>> primes(10) [2, 3, 5, 7] """ if n==2: return [2] elif n<2: return [] s=range(3,n+1,2) mroot = n ** 0.5 half=(n+1)/2-1 i=0 m=3 while m <= mroot: if s[i]: j=(m*m-3)/2 s[j]=0 while j