Figure 5.4:

Reaction-coordinate diagram.

Code for Figure 5.4

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

initdata ();

en1 = 3;
en2 = 5;

[s, V, d3trace, num] = inittrace (en1, en2);

%----------------------------
				%the following part is different from
				%morse.m, hfc.m and hfs.m
hfptable = [s, V];
save -ascii hfp.dat hfptable;
				%End of different part
%----------------------------

inittrace.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function [s, V, d3trace, num] = inittrace (en1, en2)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  R1sar = 1.29781;
  R2sar = 0.80363;
  numr = 30;

  R1r = linspace (R1sar, en1, numr)';

  opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, R2r] = ode15s (@rder, R1r, R2sar, opts);

  R1sal = 1.297;
  R2sal = 0.80363;
  numl = 1000;

  R2l = linspace (R2sal, en2, numl)';

  opts  =  odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
  [tsolver, R1l] = ode15s (@lder, R2l, R1sal, opts);

  R1rn = zeros (numr, 1);
  R2rn = zeros (numr, 1);
  num = numl+numr;

  for i = 1:numr
    R1rn(i) = R1r(numr+1-i);
    R2rn(i) = R2r(numr+1-i);
  end

  R1 = [R1rn; R1l];
  R2 = [R2rn; R2l];
  trtable = [R1, R2];

  s = zeros (num, 1);
  V = zeros (num, 1);
  D = zeros (num, 1);

  s(1) = 0;
  V(1) = Vau (R1(1), R2(1));
  D(1) = 1;

  for i = 2:num
    D(i) = i;
    s(i) = sqrt ((R1(i)-R1(i-1))^2 + (R2(i)-R2(i-1))^2) + s(i-1);
    V(i) = Vau (R1(i), R2(i));
  end

  pot = [s, V, R1, R2, D];

  for i = 1:num
    j = i*3-2;
    d3trace(j) = R1(i);
    d3trace(j+1) = R2(i);
    d3trace(j+2) = V(i);
  end

  d3trace = d3trace';

  % this seems to be unused.
  % for a = 1:3
  %   for i = 1:num
  %    j = i*3-2;
  %    d3sptrace(a,j) = R1(i)+(a-2)*0.1;
  %    d3sptrace(a,j+1) = R2(i);
  %    d3sptrace(a,j+2) = -100+(a-2);
  %   end
  % end

end%function

initdata.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function initdata ()

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  D1 = [129.3; 94.8; 129.3];
  B1 = [2.316; 2.047; 2.316];
  r10 = [0.931; 0.753; 0.931];
  D3 = [51.83; 34.56; 51.83];
  B3 = [2.229; 1.522; 2.229];
  r30 = [0.931; 0.753; 0.931];
  R_1 = [1.217; 0.443; 1.217];
  C = [1174.9; 388.8; 1174.9];
  Si = [3.090; 1.929; 3.090];
  A = [1.315; 0.762; 1.315];

end%function

Vau.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = Vau(r1,r2)
global D1 D3 B1 B3 r10 r30 R_1 C Si A
r3=r1+r2;
r(1)=r1;
r(2)=r2;
r(3)=r3;
for x = 1:3
 E1(x)=D1(x)*(exp(-2*B1(x)*(r(x)-r10(x)))-2*exp(-B1(x)*(r(x)-r10(x))));
 E3(x,1)=D3(x)*(exp(-2*B3(x)*(r(x)-r30(x)))+2*exp(-B3(x)*(r(x)-r30(x))));
 E3(x,2)=C(x)*(r(x)+A(x))*exp(-Si(x)*r(x));
  if(r(x) <= R_1(x))
   y=1;
   else
   y=2;
  end
 Q(x)=(E1(x)+E3(x,y))/2;
 J(x)=(E1(x)-E3(x,y))/2;
end
retval=Q(1)+Q(2)+Q(3)-sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))));
end%function

lder.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = lder (r2, r1)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  r3 = r1+r2;
  r(1) = r1;
  r(2) = r2;
  r(3) = r3;

  for x = 1:3
   E1(x) = D1(x)*(exp(-2*B1(x)*(r(x)-r10(x)))-2*exp(-B1(x)*(r(x)-r10(x))));
   E3(x,1) = D3(x)*(exp(-2*B3(x)*(r(x)-r30(x)))+2*exp(-B3(x)*(r(x)-r30(x))));
   E3(x,2) = C(x)*(r(x)+A(x))*exp(-Si(x)*r(x));
    if(r(x) <= R_1(x))
     y = 1;
     else
     y = 2;
    end
   Q(x) = (E1(x)+E3(x,y))/2;
   J(x) = (E1(x)-E3(x,y))/2;
  end

  %dQ(a,b) means dQ(a)/dr(b)
  for i = 1:3
     dE1(i,i) = -2*B1(i)*D1(i)*(exp(-2*B1(i)*(r(i)-r10(i)))-exp(-B1(i)*(r(i)-r10(i))));
   if(r(i) <= R_1(i))
     dE3(i,i) = -2*B3(i)*D3(i)*(exp(-2*B3(i)*(r(i)-r30(i)))+exp(-B3(i)*(r(i)-r30(i))));
   else
     dE3(i,i) = C(i)*exp(-Si(i)*r(i))*(1-Si(i)*(r(i)+A(i)));
   end
  end
  dE1(1,2) = 0;
  dE3(1,2) = 0;
  dE1(2,1) = 0;
  dE3(2,1) = 0;
  dE1(3,1) = dE1(3,3);
  dE3(3,1) = dE3(3,3);
  dE1(3,2) = dE1(3,3);
  dE3(3,2) = dE3(3,3);
  dE1(1,3) = 1000;
  dE3(1,3) = 1000;
  dE1(2,3) = 1000;
  dE3(2,3) = 1000;
  dE1(3,3) = 1000;
  dE3(3,3) = 1000;
  for i = 1:3
   for j = 1:3
    dQ(i,j) = (dE1(i,j)+dE3(i,j))/2;
    dJ(i,j) = (dE1(i,j)-dE3(i,j))/2;
   end
  end

  dV(1) = dQ(1,1)+dQ(3,1)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(2*dJ(1,1)-dJ(3,1))+J(2)*(-dJ(1,1)-dJ(3,1))+J(3)*(2*dJ(3,1)-dJ(1,1)));

  dV(2) = dQ(2,2)+dQ(3,2)-0.5/sqrt(0.5*((J(1)-J(2))*(J(1)-J(2))+(J(2)-J(3))*(J(2)-J(3))+(J(1)-J(3))*(J(1)-J(3))))*(J(1)*(-dJ(2,2)-dJ(3,2))+J(2)*(2*dJ(2,2)-dJ(3,2))+J(3)*(2*dJ(3,2)-dJ(2,2)));

  retval = dV(1)/dV(2);
end%function

rder.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

function retval = rder (r1, r2)

  global D1 D3 B1 B3 r10 r30 R_1 C Si A

  retval = 1 / lder (r2, r1, D1, D3, B1, B3, r10, r30, R_1, C, Si, A);

end%function