############################################################################## # Simulation of Nuclear Magnetic Resonance # # Reinhard A. Schumacher (All rights reserved) # Carnegie Mellon University # January 2002 # # Modifications: # 3-03-2007 - add run/stop toggle, restructure display # 3-16-2007 - improve the modeling of T1 relaxation # from visual import * from random import random, randrange from visual.graph import * print """ ############################################################################## Simulation of Nuclear Magnetic Resonance (NMR) using a classical vector model. ------------------------------------------------------------------------------ A set of spins precesses in a strong static external magnetic field. A perp- endicular oscillating magnetic field causes the reorientation of the spins, inducing a non-zero net magnetization in the plane perpendicular to the static external field. This net oscillating magnetization produces the NMR signal. Quantum physics shows that the expectation value of the quantum spin operator evolves in time in the same way as a classical spin vector. Thus, this simulation is a faithful representation of what the quantum system does. The simulation includes the effects (controllable) of a non-uniform magnetic field, the spin-spin interactions of the protons (T2), and the relaxation of the spins back to the external field direction after they have been rotated away from that axis (T1). Controls: Left-clicking spheres on the main display screen activates given action. Left-clicking anywhere else pauses/resumes animation. Left-clicking on the "Net Magnetization" panel clears the arrow "trail". Adjustments: Inside the program, various parameters can be set: numspin = number of randomly placed, randomly oriented (in phi), spins bfield0 = vector components of the strong external field dbdx = vector gradient components along X for external field dbdy = vector gradient components along Z(up!)of external field dbdz = vector gradient components along Y of external field bfield1mag = magnitude of the rotating field bfield2mag = maximum value of "random" spin-spin field: the T2 process freq1 = oscillating field frequency teeone = T1 relaxation time towards vertical axis speed = overall speed of the simulation, set depending on CPU speed R. A. Schumacher Carnegie Mellon University Version 5.0 ########################################################################### """ ########################################################################### # USER ADJUSTABLE INPUTS GO HERE: numspin = 400 #100 number of randomly placed spins bfield0 = vector(0,5,0) #0,5,0 components of the strong external field dbdx = vector(0,0.07,0) #.07 gradient components along X of external field dbdy = vector(0,0.00,0) #.05 gradient components along Z(up!)of external field dbdz = vector(0,0.00,0) #.03 gradient components along Y of external field bfield1mag = 2.0 #2. magnitude of the rotating field bfield2mag = 1.0 #0.5 maximum value of "random" spin-spin field --> T2 freq1 = 10. #10. oscillating field frequency teeone = 10. #10. T1 relaxation time towards vertical axis speed = 20 #maximum update speed per second: depends on CPU power ########################################################################### scene1 = display( title = "Nuclear Magnetic Resonance as an Ensemble of Classical Spins" , width=700, height=700, autoscale=0, background=(.3,.4,.5), x=0, y=0, forward=(-0.25,-0.25,-1)) scene2 = display(title = "NMR: Net Magnetization" , width=400, height=400, center=(0,0,0), autoscale=0, background=(.3,.4,.5), x=700, y=0, forward=(0,-.1,-1.), range=(4,4,4)) gdisplay(title='NMR: Graph', xtitle='Time', ytitle='XY Magnetization', background=(.2,.2,.2), xmax=100.0, ymax=1.1, ymin=0.0, x=700, y=400, width=500, height=300) mt = gcurve(color=color.magenta) def wirecube (s,x,y,z): c=curve (color=color.white, radius=.1) pts = [(-s+x, -s+y, -s+z),(-s+x, -s+y, s+z), (-s+x, s+y, s+z), (-s+x, s+y, -s+z), (-s+x, -s+y, -s+z), (s+x, -s+y, -s+z), (s+x, s+y, -s+z), (-s+x, s+y, -s+z), (s+x, s+y, -s+z), (s+x, s+y, s+z), (-s+x, s+y, s+z), (s+x, s+y, s+z), (s+x, -s+y, s+z), (-s+x, -s+y, s+z), (s+x, -s+y, s+z),(s+x, -s+y, -s+z)] for pt in pts: c.append(pos=pt) # # Set up the visual scene # scene1.select() #makes first scene the current default for new objects wirecube(5,0,0,0) scene2.select() wirecube(2,0,0,0) scene1.select() base = box(pos=vector(0,-5,0),height=0.1, length=10.0, width=10.0, color=(0.0,0.5,1.0)) s_reset=sphere(pos=(-5,5,5), radius=0.25, color=(0.6, 1.0, 0.6)) s_reset_label=label(pos=s_reset.pos, text="Reset Spins", yoffset=+10, height=12,opacity=0, box=0, line=0) s_osc=sphere(pos=(5,5,5), radius=0.2, color=(1.0, 0.0, 0.0)) s_osc_label=label(pos=s_osc.pos, text="Oscillating Field ON/OFF", yoffset=+10, height=12,opacity=0, box=0, line=0) s_dbdr=sphere(pos=(-5,-5,5), radius=0.25, color=(1.0, 0.0, 0.0)) s_dbdr_label=label(pos=s_dbdr.pos, text="Gradients ON/OFF", yoffset=-10, height=12,opacity=0, box=0, line=0) s_freq1dn=sphere(pos=(-5,0,5), radius=0.25, color=(0.6, 1.0, 0.6)) s_freq1dn_label=label(pos=s_freq1dn.pos, text="Lower Freq.", xoffset=-10, height=12,opacity=0, box=0, line=0) s_freq1up=sphere(pos=(5,0,5), radius=0.25, color=(0.6, 1.0, 0.6)) s_freq1dn_label=label(pos=s_freq1up.pos, text="Raise Freq.", xoffset=10, height=12,opacity=0, box=0, line=0) s_t1onoff=sphere(pos=(5,-5,5), radius=0.25, color=(1.0, 0.0, 0.0)) s_t1onoff_label=label(pos=s_t1onoff.pos, text="T1 Relax. ON/OFF", yoffset=-10, height=12,opacity=0, box=0, line=0) s_t2onoff=sphere(pos=(0,-5,5), radius=0.25, color=(1.0, 0.0, 0.0)) s_t2onoff_label=label(pos=s_t2onoff.pos, text="T2 Relax. ON/OFF", yoffset=-10, height=12,opacity=0, box=0, line=0) Spins = [] dera = pi/180. spinmag = 1.0 #magnitude of all spin vectors polarangle = 2 #a hack... for i in range(numspin): xpos = randrange(-500.,500.)/100. ypos = randrange(-500.,500.)/100. zpos = randrange(-500.,500.)/100. phiangle = randrange(0.,360.) if i==1: xpos=0. ypos=0. zpos=0. phiangle = 0. # (VPython forces us to accept "up" as the Y direction) xval = spinmag*sin(polarangle*dera)*cos(phiangle*dera) zval = spinmag*sin(polarangle*dera)*sin(phiangle*dera) yval = spinmag*cos(polarangle*dera) spinvector = vector (xval,yval,zval) posvector = vector (xpos,ypos,zpos) Spins = Spins+[arrow(pos=posvector, axis=spinvector, shaftwidth=0.10, color=color.green)] # # Draw permanent arrows representing magnetic field in space, # including the gradients. At present, we only represent grad(B(x)) non-zero. # bfield = bfield0 + (-5)*dbdx btail = -bfield.mag / 2. bleft = arrow(pos=(-5,btail,0), axis=bfield, shaftwidth=0.05, color=color.blue) bfield = bfield0 + (5)*dbdx btail = -bfield.mag / 2. brght = arrow(pos=( 5,btail,0), axis=bfield, shaftwidth=0.05, color=color.blue) btail = -bfield0.mag / 2. bback = arrow(pos=( 0,btail,-5),axis=bfield0, shaftwidth=0.05, color=color.blue) bfrnt = arrow(pos=( 0,btail,5), axis=bfield0, shaftwidth=0.05, color=color.blue) # brot = arrow(pos=( 0,0,0), axis=(bfield1mag,0,0), shaftwidth=0.05, color=color.cyan, visible=0,display=scene1) bosc = arrow(pos=( 0,0,0), axis=(bfield1mag,0,0), shaftwidth=0.05, color=color.magenta, visible=0,display=scene1) bosc2 = arrow(pos=( 0,0,0), axis=(bfield1mag,0,0), shaftwidth=0.05, color=color.magenta, display=scene2, visible=0) netmag = arrow(pos=(0,0,0), axis=(0,0,0), shaftwidth=0.2, color=color.green, display=scene2) trail = curve(display=scene2) #follows the net magentization in scene2 btrail = curve() #follows the rotating field vector tmark = sphere(pos = (-5,-5,5.2), radius = 0.05, color = color.red, visible=0 ) # # Set up the physics # tflag = 0 t1fact = 0 #start with no relaxation towards B_zero t2fact = 0 #start with no spin-spin interaction gradfact = 0 #start with external field gradients turned off gamma = 1.0 magmoment = 2.0 ratio = magmoment/gamma freq0 = ratio*bfield0.mag #Larmor frequency at center of space t = 0. #global time t1 = 0. #time since rotating field was turned on dt = 0.03 #time step print "Larmor frequency=",freq0,"cyles/time, Rotating field frequency=",freq1,"cycles/time" while 1: rate(speed) #VPython's maximum update rate in inverse seconds t=t+dt if tflag==1: #this flag says the rotating field is "ON" t1 = t1 + dt angle1 = freq1*t1 #angle of rotating B-field in horizontal plane brot.axis.x = bfield1mag*cos(angle1) #effective rotating field brot.axis.z = -bfield1mag*sin(angle1) bosc.axis.x = brot.axis.x #actual applied oscillating field... bosc2.axis.x= brot.axis.x #...duplicate for 2nd scene bfield = bfield0 + 0.5*bosc.axis if angle1<2.5*pi: btrail.append(pos=brot.axis) elif tflag==0: bfield = bfield0 sumx = 0 sumy = 0 sumz = 0 for i in range(numspin): xxx = Spins[i].pos.x yyy = Spins[i].pos.y zzz = Spins[i].pos.z btot = bfield + xxx*dbdx*gradfact #increment field according to gradients btot = btot + yyy*dbdy*gradfact #increment field according to gradients btot = btot + zzz*dbdz*gradfact #increment field according to gradients if t2fact==1: branx = 0 # bfield2mag*(2.*random()-1.) brany = bfield2mag*(2.*random()-1.) #Let the T2 process act only in XY branz = 0 # bfield2mag*(2.*random()-1.) brandom = vector(branx,brany,branz) #simulates spin-spin "T2" processes btot = btot + brandom #add a different random field at each spin danglerad = ratio*(btot.mag)*dt #increment in precession angle Spins[i].rotate(angle=danglerad, axis=btot) # # Simulate either a smooth exponential T1 decay back to original direction, # or a stochastic step-wise decays process. Either way, we assume the # process is an unspecified non-radiative small-step process that incrementally # changes the expectation value of the spin-up/spin-down ratio of each spin # separately. Assume that the T1 process has no effect on the XY phase of the # spin wavefunctions. # sxaxis = Spins[i].axis.x #let the spins decay back to external syaxis = Spins[i].axis.y #field direction with "T1" constant szaxis = Spins[i].axis.z # sdy = +(1.-syaxis)*dt/teeone * t1fact #smooth if random() < dt/teeone*t1fact: #stochastic sdy = +0.4 #arbitrary selection of "energy" scale for T1 interaction else: sdy = 0.0 syaxis = syaxis + sdy if syaxis>1: syaxis = 1. if syaxis<-1: syaxis = -1. sxzmax = sqrt(1. - syaxis**2.) sxz = sqrt(sxaxis**2. + szaxis**2.) if sxz > 0.: sizrat = sxzmax/sxz else: sizrat = 10. #don't scale transverse projections if sizrat < 1 & t1fact==1: sxaxis = sxaxis*sizrat szaxis = szaxis*sizrat Spins[i].axis.x = sxaxis Spins[i].axis.y = syaxis Spins[i].axis.z = szaxis redfrac = 1. - (syaxis+1.)/2. grnfrac = 1. - redfrac blufrac = 0. Spins[i].color=(redfrac,grnfrac,blufrac) sumx = sumx+sxaxis sumy = sumy+syaxis sumz = sumz+szaxis # # Make the net magnetization arrow some blended color... # blufrac = 0. redfrac = 1.-((sumy/numspin)+1.)/2. grnfrac = 1.-redfrac netmag.color =(redfrac,grnfrac,blufrac) netmag.axis =(2*sumx/numspin,2*sumy/numspin,2*sumz/numspin) xymag = sqrt(netmag.axis.x**2 + netmag.axis.z**2)/2. mt.plot(pos=(t,xymag)) trail.append(pos=netmag.axis) # # Mousing around... # if scene1.mouse.events: m = scene1.mouse.getevent() # # Click on a labeled sphere to get the stated action. # if m.click: if m.pick is s_reset: print "Resetting spins..." for i in range(numspin): Spins[i].axis.x = (2.*random()-1.)*0.01 Spins[i].axis.y = spinmag Spins[i].axis.z = (2.*random()-1.)*0.01 Spins[i].color = color.green trail.visible = 0 trail = curve(display=scene2) mt.visible = 0 mt = gcurve(color=color.cyan) btrail.visible = 0 btrail = curve() t = 0 t1 = 0 elif m.pick is s_dbdr: if gradfact == 0: gradfact = 1 s_dbdr.color = color.green print "External field gradients turned ON" else: gradfact = 0 s_dbdr.color = color.red print "External field gradients turned OFF" elif m.pick is s_freq1up: freq1 = freq1 + 2. print "Raising oscillation frequency..." print "Larmor frequency=",freq0,"cyles/time, Rotating field frequency=",freq1,"cycles/time" elif m.pick is s_freq1dn: freq1 = freq1 - 2. print "Lowering oscillation frequency..." print "Larmor frequency=",freq0,"cyles/time, Rotating field frequency=",freq1,"cycles/time" elif m.pick is s_t1onoff: if t1fact == 0: t1fact = 1 s_t1onoff.color = color.green print "T1 relaxation is turned ON" else: t1fact = 0 s_t1onoff.color = color.red print "T1 relaxation is turned OFF" elif m.pick is s_t2onoff: if t2fact == 0: t2fact = 1 s_t2onoff.color = color.green print "T2 relaxation is turned ON" else: t2fact = 0 s_t2onoff.color = color.red print "T2 relaxation is turned OFF" elif m.pick is s_osc: if tflag ==0: tmark.visible = 1 btrail = curve() btrail.visible = 1 bosc.visible = 1 bosc2.visible = 1 brot.visible = 1 tflag = 1 t1 = 0. #start the rotating field at phi=0 s_osc.color = color.green print "Rotating field turned ON" scene1.title = 'NMR: Rotating field turned ON' else: tmark.visible = 0 btrail.visible = 0 bosc.visible = 0 bosc2.visible = 0 brot.visible = 0 tflag = 0 s_osc.color = color.red print "Rotating field turned OFF" scene1.title = 'NMR: Rotating field turned OFF' else: print "Pausing animation" scene1.mouse.getclick() #inhibits animation until the mouse is clicked print "Resuming animation" if scene2.mouse.events: m = scene2.mouse.getevent() print "Clearing trail" trail.visible = 0 trail = curve(display=scene2) #follows the net magentization in scene2 # else: # print "Pausing animation" # scene2.mouse.getclick() #inhibits the calculation until the mouse is clicked # print "Resuming animation" # END OF CODE