Chem549 Assignment 6
The Split Operator Method with time-dependent potentials.
For the anharmonic oscillator (particle in a Morse potential):
clearall;
L=1024;
h=6.62606957e-34;
hbar=h/2/pi;
dt=0.025e-19;
m=9.10938e-31;
dx=1e-12;
i=sqrt(-1);
length=L*dx;
ke=7000;
De=1e-16;
a=sqrt(ke/2/De);
lam=sqrt(2*m*De)/a/hbar;
v0=a/2/pi*sqrt(2*De/m);
for j=1:L
x(j)=dx*(j);
psi_x(j)=exp(-(x(j)-.22e-9)^2/2e-21);
p(j)=2*(j-1)*pi/length;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
v(j)=De*((1-exp(-a*(x(j)-.15e-9)))^2-1)+De;
opk(j)=exp(-i*dt*hbar*hbar*p(j)^2/4/hbar/m);
opv(j)=exp(-i*dt*v(j)/hbar);
end;
norm=sum(abs(psi_x));
psi_x=psi_x/norm;
psi_initial=psi_x;
n=0;
eigen1=0;
eigen2=0;
eigen3=0;
eigen4=0;
for z=1:1600000
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
psi_x=opv.*psi_x;
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
p(z)=sum(conj(psi_x).*psi_initial);
eigen1=eigen1+exp(i*4.638248699000001e-18*dt*z/hbar)*psi_x;
eigen2=eigen2+exp(i*1.341779087925000e-17*dt*z/hbar)*psi_x;
eigen3=eigen3+exp(i*2.173350818960000e-17*dt*z/hbar)*psi_x;
eigen4=eigen4+exp(i*2.981731306500000e-17*dt*z/hbar)*psi_x;
if (mod(z,100)==0)
n=n+1;
remember(:,n)=psi_x;
end;
end;
eigen1=eigen1/sum(abs(eigen1));
eigen2=eigen2/sum(abs(eigen2));
eigen3=eigen3/sum(abs(eigen3));
eigen4=eigen4/sum(abs(eigen4));
I know there are a lot of parameters in the beginning; these are just to correctly code the Morse potential which is a bit complicated.
The program above is designed to calculate the lowest eigenstates of the harmonic potential; also note that I get the energy of the state first and then plug that into:
eigen1=eigen1+exp(i*4.638248699000001e-18*dt*z/hbar)*psi_x;
and then re-run the calculation again. For convenience, I have already included the correct energy.
I can also calculate the energies of the next few states:
1.341779087925000e-17
2.173350818960000e-17
2.981731306500000e-17
And I will save the first eigenstate to the hard drive:
> save eigen1 eigen1
Here is what the first two look like in the Morse potential:
Next, we want to use theground state as the initial wavefunction, so I altered the code as follows:
clearall;
L=1024;
h=6.62606957e-34;
hbar=h/2/pi;
dt=0.075e-19;
m=9.10938e-31;
dx=1e-12;
i=sqrt(-1);
length=L*dx;
ke=7000;
De=1e-16;
a=sqrt(ke/2/De);
lam=sqrt(2*m*De)/a/hbar;
v0=a/2/pi*sqrt(2*De/m);
for j=1:L
x(j)=dx*(j);
p(j)=2*(j-1)*pi/length;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
v(j)=De*((1-exp(-a*(x(j)-.15e-9)))^2-1)+De;
opk(j)=exp(-i*dt*hbar*hbar*p(j)^2/4/hbar/m);
opv(j)=exp(-i*dt*v(j)/hbar);
end;
loadeigen1;
psi_x=eigen1;
norm=sum(abs(psi_x));
psi_x=psi_x/norm;
n=0;
for z=1:133150
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
psi_x=opv.*psi_x;
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
if (mod(z,50)==0)
n=n+1;
remember(:,n)=psi_x;
end;
end;
Note the differences- I don’t make up some initial wavefunction, I read it in which is why I saved it before. So go ahead and cut and paste and then do the following to visualize:
for v=1:n plot(x,abs(remember(:,v))); axis([-.5e-10 5e-10 0 .015]); M(v)=getframe; end;
Note how the eigenstatehardly changes over time. What slight fluctuation you may see is due to the inherent error in the method.
Now if I’m going to simulate a time varying electric field (a photon) simply note that I’m doing the same thing as putting this quantum entity in a slanty shanty that pivots like a see-saw. Here is how I will do so:
clearall;
L=1024;
h=6.62606957e-34;
hbar=h/2/pi;
dt=0.075e-19;
m=9.10938e-31;
dx=1e-12;
i=sqrt(-1);
length=L*dx;
ke=7000;
De=1e-16;
a=sqrt(ke/2/De);
lam=sqrt(2*m*De)/a/hbar;
v0=a/2/pi*sqrt(2*De/m);
for j=1:L
x(j)=dx*(j);
p(j)=2*(j-1)*pi/length;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
opk(j)=exp(-i*dt*hbar*hbar*p(j)^2/4/hbar/m);
end;
loadeigen1;
psi_x=eigen1;
norm=sum(abs(psi_x));
psi_x=psi_x/norm;
psi_initial=psi_x;
n=0;
freq=7.547169811320758e-17/dt;
for z=1:133150
v(1:L)=De*((1-exp(-a*(x(1:L)-.15e-9))).^2-1)+De+1.25e-8.*(x(1:L)-.15e-9)*sin(2*pi*z/freq);
opv(1:L)=exp(-i*dt*v(1:L)/hbar);
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
psi_x=opv.*psi_x;
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
if (mod(z,50)==0)
n=n+1;
remember(:,n)=psi_x;
end;
end;
First note that there was no need to define the potential initially as it (and the operator opv) are now in the time loop section. Note that I have the same old Morse potential, but I added this term:
1.25e-8.*(x(1:L)-.15e-9)*sin(2*pi*z/freq);
The first number is an approximate representation of the strength of a photon field hitting the particle in the oscillator. The second part is the “slant” centered at the middle of the potential. The third term, sine, causes the potential to reverse directions according to the frequency. The frequency is derived from the energies of the states as follows
1.341779087925000e-17-4.638248699000001e-18
ans =
8.779542180249999e-18
ans/h
ans =
1.325000000000000e+16
Now that is the frequency that corresponds to the energy difference between the ground and first excited states. To translate it to time, I calculate 1/frequency:
> 1/ans
ans =
7.547169811320756e-17
Note this is where the first term in the expansion comes from. Now I divide by dt which represents the number of time steps necessary to complete a full cycle (as dictated by the time loop variable z).
Now run the program and let’s visualize it.
> for v=1:n plot(x,abs(remember(:,v))); axis([-.5e-10 5e-10 0 .015]); M(v)=getframe; end;
In the beginning, you see the ground state:
But by the end of the trajectory, you have the first excited state (or something close to it):
Sometimes you see a third peek. This is due to the fact that the transition energy (same as frequency) from the 1st excited state to the 2nd isn’t very different (it’s a little bit lower) so some of the 2nd excited state character is making an appearance.
Now run the trajectory twice as long, and the last thing you see is this:
Why are we back to the ground state? Because you just discovered stimulated emission. Yes, this program proves lasers work.
Last bit, plug in the energy difference between the ground state and 2nd excited state like I did above. Calculate the frequency, and drive the system to make that transition. Doesn’t do much does it? After all when you doubled the z loop you saw the wavefunction go from the 1st excited state and back to ground (two transitions). You don’t really observe even one transition. Is there something wrong? Try tripling the photon field strength to 10e-8, and the transition from ground, to 2nd excited and back again is observed. Congrats- you just discovered the vibrational overtone. Now you have to figure out why it takes so long.
The actual problem set.
For this one, we will work with the particle in a box model from last time. Go back and make some changes to the program for the particle-in-a-box question, and then calculate the first eigenstate (ground state).
clearall;
L=1024;
hbar=1.0545717e-34;
dt=0.05e-19;
m=9.10938e-31;
dx=1e-12;
i=sqrt(-1);
length=L*dx;
height=6e-6;
for j=1:L
x(j)=dx*(j-L/2);
psi_x(j)=exp(-(x(j)+0.3e-10)^2/2e-21);
p(j)=2*(j-1)*pi/length;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
v(j)=height;
if (j>L/3&j<2*L/3)
v(j)=0;
end;
opk(j)=exp(-i*dt*hbar*hbar*p(j)^2/4/hbar/m);
opv(j)=exp(-i*dt*v(j)/hbar);
end;
norm=sum(abs(psi_x));
psi_x=psi_x/norm;
psi_initial=psi_x;
n=0;
eigen1=0;
eigen2=0;
eigen3=0;
eigen4=0;
for z=1:1600000
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
psi_x=opv.*psi_x;
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
p(z)=sum(conj(psi_x).*psi_initial);
eigen1=eigen1+exp(i*4.969552058105549e-19*dt*z/hbar)*psi_x;
eigen2=eigen2+exp(i*2.070646690877312e-18*dt*z/hbar)*psi_x;
eigen3=eigen3+exp(i*4.638248587565180e-18*dt*z/hbar)*psi_x;
eigen4=eigen4+exp(i*8.282586763509248e-18*dt*z/hbar)*psi_x;
if (mod(z,100)==0)
n=n+1;
remember(:,n)=psi_x;
end;
end;
eigen1=eigen1/sum(abs(eigen1));
eigen2=eigen2/sum(abs(eigen2));
eigen3=eigen3/sum(abs(eigen3));
eigen4=eigen4/sum(abs(eigen4));
end;
Save the ground state to the computer.
1.Now here is some code that can perturb the particle in a box. Calculate the first three transition timescales (see where I left a question mark?):
clearall;
L=1024;
h=6.62606957e-34;
hbar=h/2/pi;
dt=0.05e-19;
m=9.10938e-31;
dx=1e-12;
i=sqrt(-1);
length=L*dx;
height=6e-6;
for j=1:L
x(j)=dx*(j-L/2);
p(j)=2*(j-1)*pi/length;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
if (j>L/2)
p(j)=2*(L-j+1)*pi/length;
end;
opk(j)=exp(-i*dt*hbar*hbar*p(j)^2/4/hbar/m);
end;
v(1:L)=height;
loadeigen1;
norm=sum(abs(eigen1));
psi_x=eigen1/norm;
n=0;
%freq1=?/dt;
%freq2=?/dt;
%freq3=?/dt;
v(1:L)=height;
opv(1:L)=exp(-i*dt*v(1:L)/hbar);
for z=1:300000
v(342:682)=8e-9*x(342:682)*sin(2*pi*z/freq1);
opv(342:682)=exp(-i*dt*v(342:682)/hbar);
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
psi_x=opv.*psi_x;
psi_k=fft(psi_x);
psi_k=opk.*psi_k;
psi_x=ifft(psi_k);
psi_x=psi_x/sum(abs(psi_x));
p(z)=sum(conj(psi_x).*psi_initial);
if (mod(z,400)==0)
n=n+1;
remember(:,n)=psi_x;
end;
end;
Send me the frequencies. Drive the system at the first frequency and tell me what happens. Then the second frequency- it’s very different. Next, increase the length of the trajectory by a factor of 8:
for z=1:300000*8 …
and drive the system with the 3rd frequency. Similar to the first, but took much longer. Note you want to visualize the results with the following:
for v=1:n plot(x,abs(remember(:,v))); axis([-5e-10 5e-10 0 .0075]); M(v)=getframe; end;
2.Transition intensity can be calculated as follows:
Since you know the first four eigenstates of the particle in the box, use the symbolic integrator from Mathematica to calculate the above where i is the ground state the f are the first three excited states. Negative numbers don’t matter- just their absolute magnitude. Now do you see why nothing happened when you drove the system at the 2nd frequency, and why the third frequency required much longer to make the transition?