Mathematica Noteook for Calculating Rate of Excitation Transfer between two BChls

NOTE: To execute a group of command lines, click anywhere within a cell covered by a blue bracket on the right which has a semi-arrow, press and hold the Shift key and then press Enter.

STEP 1: Define constants used in this calculation and Load a package called "VectorAnalysis".

const :=  505644;   
hbar := 6.582*(10^(-16))*8065.82*(10^12);    
J := 10.5*(10^(-4));    
J1 :=  2.785*(10^(-4));
J2 := 6.295*(10^(-10));
<< Calculus`VectorAnalysis`

STEP 2: Compute the directional vector d1 of the transition dipole moment  of the BChl B850a in segment BCA1.

rNB1  =  {xNB1,  yNB1,  zNB1};
rND1  =  {xND1,  yND1,  zND1};
D1  =  rND1 - rNB1;
D1sq = DotProduct[D1, D1];
d1 = D1/Sqrt[D1sq]

STEP 3: Compute the directional vector d3 of the transition dipole moment of the BChl B850a in segment BCA2.

rNB3 = {xNB3,  yNB3,  zNB3};
rND3 = {XND3,  yND3,  zND3};
D3  =  rND3 - rNB3;
D3sq = DotProduct[D3, D3];
d3 = D3/Sqrt[D3sq]

STEP 4: Compute the separation distance L13 between the two central Mg2+ ions of the two BChls B850a in segments BCA1 and BCA2.

rMG1 = {xMG1,  yMG1,  zMG1};
rMG3 = {xMG3,  yMG3,  zMG3};
R13 = rMG1 - rMG3;
R13sq = DotProduct[R13, R13];
L13 = Sqrt[R13sq]

STEP 5: Calculate the directional vector of the relative position vector r13 between the two Mg2+ ions

r13 = R13/L13

STEP 6: Calculate the induced dipole-induced dipole interaction term U13

d1d3 = DotProduct[d1, d3];
d1r13 = DotProduct[d1, r13];
d3r13 = DotProduct[d3, r13];
U13 = 
      const*(d1d3 - 3*d1r13*d3r13)/(L13^3)

STEP 7: Calculate the excitation transfer rate k13 between two nearest BChls B850a

k13 = ((2*
              Pi)/hbar)*(U13^2)*J

STEP 8: Calculate the forward rate kf of excitation transfer from BChl B800 (Segment BCA1)  to BChl B850a (segment BCA1)

sNB1  =  {27.736, -11.012,  
          89.716};
sND1  =  {24.434, -13.149,  
          90.602};
Ds1  =  sND1 - sNB1;
Ds1sq = DotProduct[Ds1, Ds1];
ds1 = Ds1/Sqrt[Ds1sq];
sMG1 = {26.069, -12.126,  90.141};

    
Rs11 = sMG1 - rMG1;
Rs11sq = DotProduct[Rs11, Rs11];
Ls11 = Sqrt[Rs11sq];

rs11 = Rs11/Ls11;
ds1d1 = DotProduct[ds1, d1];
ds1rs11 = DotProduct[ds1, rs11];
d1rs11 = DotProduct[d1, rs11];
Us11 = 
        const*(ds1d1 - 3*ds1rs11*d1rs11)/
            (Ls11^3);
kf = ((2*
              Pi)/hbar)*(Us11^2)*J1

STEP 9: Compute the backward rate kb of excitation transfer from BChl B850a (segment BCA1)  to B800 (segment BCA1)

kb = ((2*
              Pi)/hbar)*(Us11^2)*J2