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<half:
s[j]=0
j+=m
i=i+1
m=2*i+3
return [2]+[x for x in s if x]
number_to_try = 1000000
list_of_primes = primes(number_to_try)
print list_of_primes[10001]
cProfile.run('primes(1000000)')