function ret = FMO
% jacob.biamonte@comlab.ox.ac.uk 
% should produce fig 4 from:
% A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U.S.A. 106, 17255 (2009)

hbar = 1.055e-34;
Ham =  ...
    [12410, -87.7, 5.5, -5.9, 6.7, -13.7, -9.9; ...
    -87.7, 12530, 30.8, 8.2, 0.7, 11.8, 4.3; ...
    5.5, 30.8, 12210, -53.5, -2.2, -9.6, 6.0; ...
    -5.9, 8.2, -53.5, 12320, -70.7, -17.0, -63.3; ...
    6.7, 0.7, -2.2, -70.7, 12480, 81.1, -1.3; ...
    -13.7, 11.8, -9.6, -17.0, 81.1, 12630, 39.7; ...
    -9.9, 4.3, 6.0, -63.3, -1.3, 39.7, 12440];

V=diag(diag(Ham));

K = V-Ham;

 [X, lambda] = eig (Ham);
 E=diag(lambda)
 %X=orth (X, .00001);

% total time, in the paper they use fs, 1e-15
Total = 1000*1e-15;
T = 0:Total/1000:Total;

p1 = zeros(1,length(T));
p2 = zeros(1,length(T));
p3 = zeros(1,length(T));
p4 = zeros(1,length(T));
p5 = zeros(1,length(T));
p6 = zeros(1,length(T));
p7 = zeros(1,length(T));

one= [1, 0, 0, 0, 0, 0, 0]';
two= [0, 1, 0, 0, 0, 0, 0]';
three= [0, 0, 1, 0, 0, 0, 0]';
four= [0, 0, 0, 1, 0, 0, 0]';
five= [0, 0, 0, 0, 1, 0, 0]';
six= [0, 0, 0, 0, 0, 1, 0]';
seven= [0, 0, 0, 0, 0, 0, 1]';

psi = one;%+seven+three+four;%+two+three+four+five+six+seven; 
psi = psi*1/norm(psi);

% NOTE: need to convert everything to the same units, so that the argument of expm is dimensionless.
% NOTE: One way is to convert time to cm  that  is inverse wavenumbers. Have a look at the following Mathematica function:
% in mathematica:  TimeInverseWavenumbers[seconds_] := N[2*Pi*(3*10^10)*seconds];


for w = 1:length(T)
        psi = expm(-i*(K)*T(w)*1.98645e8)*psi; %here the extra numbers convert time to cm so expm is dimensionless  
        psi = psi*1/norm(psi); %should do nothing unless something bad happens 
        p1(w) = abs(dot (one,psi))^2;
        p2(w) = abs(dot (two,psi))^2;
        p3(w) = abs(dot (three,psi))^2;
        p4(w) = abs(dot (four,psi))^2;
        p5(w) = abs(dot (five,psi))^2;
        p6(w) = abs(dot (six,psi))^2;
        p7(w) = abs(dot (seven,psi))^2;
    end;

close all; 
figure;
hold on;
plot(p1','--r','LineWidth',1.5)
plot(p2','--g','LineWidth',1.5);
plot(p3','--b','LineWidth',1.5);
plot(p4','--m','LineWidth',1.5);
plot(p5','--c','LineWidth',1.5);
plot(p6','--y','LineWidth',1.5);
plot(p7','--k','LineWidth',1.5);
 
 legend('1','2','3','4','5','6','7')
 title('Time evolution of the population at each BChl of the FMO complex')
ylabel('Population of BChl')
xlabel('Time (fs)')
 
print('-dpsc','BChl1.ps')

 
 % NOTE: watch the units, as they need to be adjusted slightly! 
 
 
 
 