% written by Jim Faeder for Quantum Molecular Dynamics, Weizmann, 2000.
% This program accompanies Ch. 13 of Introduction to Quantum Mechanics:
% A Time Dependent Perspective, by David Tannor (University Science Books, 2001).

% function p_final=i2_interfere( phi, Tau)   

% Example of wavepacket interferometry in iodine molecule.  Computes the first
% order wavefunction as a function of time for two Gaussian pulses of width
% FWHM separated by time-delay Tau (in fs).  
% If plot_wave>0, the first order wavefunction (absolute value in green, real
% part in red) is plotted during the wavepacket propagation.  At the end of 
% the propagation, the electric field intensity, the excited state population, 
% and the expectation value of the bond length in the excited state are
% plotted as a function of time in the second window.  The function returns
% the final population in the excited state. 
 
% define various conversion constants
AU_TO_EV=27.2113961;
NM_TO_AU=45.56335;
AU_TO_CM=2.194746e5;
AU_TO_ANG=0.529177249;
AU_TO_FS=0.0241889;
AMU_TO_AU=1.8228885e3;
conv_r= AU_TO_ANG;
conv_en=AU_TO_CM;
conv_t= AU_TO_FS;

% Turn on and off plotting of the wavepacket wavefunction
plot_wave=1;

% set up coordinate grid
rmin= 4.3; % Note: in Bohr
rmax=  9.0; 
ngrid=128; % Should be power of 2 for good performance of FFT
dr= (rmax-rmin)/(ngrid-1);
r= rmin:dr:rmax;

% mass of I2
mass= (126.904473/2)*AMU_TO_AU;

% X (ground) state of I2---not used here for dynamics but useful for
% visualization 
De_X= 12440/AU_TO_CM;
be_X= 1.875*AU_TO_ANG;
re_X= 2.656/AU_TO_ANG;
we_X= be_X*sqrt(2*De_X/mass);

%V1= 1/2*mass*(we_X*(r-re_X).^2);
V1= De_X*(1-exp(-be_X*(r-re_X))).^2;

% B (excited) state of I2
De_B=  5169/AU_TO_CM;
be_B= 1.696*AU_TO_ANG;
re_B= 3.025/AU_TO_ANG;
Te_B= 15771/AU_TO_CM;
we_B= be_B*sqrt(2*De_B/mass);

%V2= 1/2*mass*(we_B*(r-re_B).^2) + Te_B;
V2= De_B*(1-exp(-be_B*(r-re_B))).^2 + Te_B;

% initial ground state wavefunction (ground state of harmonic oscillator---an
% approximation here, but a good one.)
w_vib= (210)/AU_TO_CM;
beta= sqrt(mass*w_vib);
psi_init= exp(-0.5*beta^2*((r-re_X).^2));
E_init= w_vib/2;
norm= 1/sqrt(dot(psi_init, psi_init));
psi_init= norm*psi_init;

% pulse parameters
FWHM=  50.0/AU_TO_FS; % This is FWHM of pulse intensity, not amplitude
Tau = 200.0/AU_TO_FS; % Time delay between pulse centers
w0_1= 17000/AU_TO_CM; % Pulse 1 frequency
w0_2= w0_1;
t0_1=  2*FWHM; % center of laser pulse 1
t0_2=  t0_1 + Tau;
phi_2= 0; % phase of second pulse relative to first pulse
% Pulse intensities in atomic units
e0_1= 0.00386;     
e0_2= 0.00386;     
% Transition dipole in atomic units
mu12= 0.181;    

% integration parameters
dt      =     0.1 /AU_TO_FS;
t_sample=    10   /AU_TO_FS;
t_end   =   t0_2 + 300/AU_TO_FS;
update_freq= 1; % update_freq*t_sample is the interval at which plots updated
n_steps_sample= round(t_sample/dt);
n_sample= round(t_end/t_sample);
t_end= n_sample*t_sample;

% pulse function
t= 0:dt:t_end;
alpha= 2*log(2)/FWHM^2;
E= e0_1*exp(-alpha*(t-t0_1).^2).*exp(-i*w0_1*(t-t0_1)) + ...
   e0_2*exp(-alpha*(t-t0_2).^2).*exp(-i*w0_2*(t-t0_2))*exp(i*phi_2);

cpu_start= cputime;
clock_start= clock;

% set up k grid
dk= 2*pi/(ngrid*dr);
k=[[0:ngrid/2-1],[-ngrid/2:-1]]*dk;

% set up kinetic energy and k-space propagator
kin= (0.5*(k.^2)/mass).';

% set up propagators
u_kin= exp(-i*dt*kin);
u_kin_half= exp(-i*0.5*dt*kin);
u_pot1= exp(-i*dt*(V1.'));
u_pot2= exp(-i*dt*(V2.'));

psi_0=    psi_init.';          % unperturbed wavefunction
psi_1=    zeros(ngrid, 1);     % First order wavefunction

% Plot potential curves and initial wavefunction
if (plot_wave)
   figure(1);
   subplot(1,1,1);
   clf;
   pot_plot=plot( r*conv_r, [V1*conv_en; V2*conv_en]);
   axis manual;
   axis([conv_r*min(r) conv_r*max(r) 0 conv_en*V2(ngrid)]);
   hold on;
   wave_scale= 0.2*De_X*conv_en/(max(psi_init)-min(psi_init));
   psi_plot=plot( r*conv_r, wave_scale*psi_init+ E_init*conv_en, ...
                  r*conv_r, 0*psi_init+ (E_init+w0_1)*conv_en,...
                  r*conv_r, 0*psi_init+ (E_init+w0_1)*conv_en);
end

tot_steps=0;
t=0;
n_sample=0;
max_steps= round(t_end/dt);
pop_1=zeros(1,t_end/t_sample)*inf;%-----------------------------------
E1=zeros(1,(t_end/dt)+1)*inf;%----------------------------------------


in = 1; %------------------------------------------
h1=figure(1);
set(h1,'position',[0, 50,400, 400])%----------------------------
while (tot_steps<max_steps)
   % propagate wavefunctions 1/2 step forward
   U= u_kin_half;
   psi_1=    ifft(U.*fft(psi_1));

   % Propagate until next sample time 
   for step=1:n_steps_sample
      % potential step
      psi_0=    psi_0*exp(-i*dt*E_init);
      psi_1=    u_pot2.*psi_1;
      % kinetic step
      if (step==n_steps_sample)
         % apply only 1/2 kinetic step when wavefunction is sampled
         U= u_kin_half;
      else
         U= u_kin;
      end
      psi_1=    ifft(U.*fft(psi_1));

      % Add contribution due to field interactions
      W1= E(tot_steps+1)*mu12;
      psi_1=    psi_1    + i*dt*(W1.*psi_0);
      t=t+dt;
      tot_steps=tot_steps+1;
       E1(tot_steps+1)=E(tot_steps+1)+.0000001;%--------------------------------
   end
  
   n_sample=n_sample+1;
   pop= dot(psi_1,psi_1);
   pop_1(n_sample)= pop;
   if (pop)
      r_1(n_sample)= real(dot(psi_1, (r.').*psi_1)/pop_1(n_sample));
   end
   if (plot_wave & mod(tot_steps/n_steps_sample, update_freq)==0)
       %s= sprintf('\\psi_1: Time= %.0f fs', conv_t*t);
       %title(s);
       set(psi_plot(2), 'YData', 10*wave_scale*abs(psi_1)+ ...
                                   (E_init+w0_1)*conv_en);
      set(psi_plot(3), 'YData', 10*wave_scale*real(psi_1)+ ...
                                   (E_init+w0_1)*conv_en);
                               
    
  
     M1(in) = getframe(h1); %----------------------------------
      
  
      
h2=figure(2);
set(h2,'position',[500, 500,400, 400])%----------------------------
clf;
subplot(2,1,1);
%s= sprintf('\\psi_1: Time= %.0f fs', conv_t*t);
%title(s);

plot((0:dt:t_end)*conv_t, abs(E1).^2);%---------------------------------------

%ylabel('Field intensity','fontsize',13);
title('Field intensity');     
axis([0 600 0 1.5*10^-5]);

subplot(2,1,2);
plot((t_sample:t_sample:t_end)*conv_t, pop_1);
%ylabel('Excited state population','fontsize',13);
title('Excited state population');      
axis([0 600 0 0.015]);

M2(in) = getframe(h2); %----------------------------------
in = in+1;  %----------------------------------------

  end
end
%movie2avi(M1,'i2_intrfr1.avi','compression','cinepak');
%movie2avi(M2,'i2_intrfr2.avi','compression','cinepak');

return;