Figure 5.3:

Contour representation of the potential-energy surface.

Code for Figure 5.3

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
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
% 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.
%
% Revised 8/22/2018

p = initdata();

an1 = 0;
en1 = 3;

an2 = 0;
en2 = 3;

n1 = 301;
n2 = 301;

step1 = (en1 - an1) / (n1 - 1);
step2 = (en2 - an2) / (n2 - 1);

d3table = zeros(n2, 3*n1);
table = zeros(n1, n2);

idx = [1, 2, 3];

for p1 = 1:n1
    r1 = (p1 - 1)*step1 + an1;

    for p2 = 1:n2
        r2 = (p2 - 1)*step2 + an2;
        tmp = Vau(r1, r2, p);
        d3table(p2,idx) = [r1, r2, tmp];
        table(p1,p2) = tmp;
    end

    idx = idx + 3;
end

x = d3table(1,1:3:end)';
y = d3table(:,2);
z = d3table(:,3:3:end);

levels = [3000, 1000, 200, -20, -40, -60, -80, ...
         -90.5, -92, -95, -100, -110, -120];

c = contourc(x, y, z, levels);

[fid, errmsg] = fopen('hfc.dat', 'w');
fprintf(fid, '% [contours]\n');

idx = 1;
nel = numel(c);

while (idx <= nel)
    fprintf(fid, '% contour level %f\n', c(idx));
    idx = idx + 1;
    len = c(idx);
    fprintf(fid, '%f %f\n', c(idx+1:idx+2*len));
    fprintf(fid, '\n');
    idx = idx + 2 * len + 1;
end

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

fprintf(fid, '\n% [min]\n');

for i = 1:num
    j = i*3-2;
    fprintf(fid, '%f %f %f\n', d3trace(j:j+2));
end

fclose(fid);

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
38
39
40
41
42
43
44
45
% 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.
%
% Revised 8/22/2018

function retval = Vau(r1, r2, p)
    r3 = r1 + r2;
    r(1) = r1;
    r(2) = r2;
    r(3) = r3;

    for x = 1:3
        E1(x) = p.D1(x)*(exp(-2*p.B1(x)*(r(x) - p.r10(x))) - ...
                    2*exp(-p.B1(x)*(r(x) - p.r10(x))));
        E3(x,1) = p.D3(x)*(exp(-2*p.B3(x)*(r(x) - p.r30(x))) + ...
                    2*exp(-p.B3(x)*(r(x) - p.r30(x))));
        E3(x,2) = p.C(x)*(r(x) + p.A(x))*exp(-p.Si(x)*r(x));

        if (r(x) <= p.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

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
% 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 p = initdata()

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

end%function

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
% 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.
%
% Revised 8/22/2018

function [s, V, d3trace, num] = inittrace(en1, en2, p)
    R1sar = 1.29781;
    R2sar = 0.80363;
    numr = 30;

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

    opts = odeset('AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
    [tsolver, R2r] = ode15s(@(r1, r2) rder(r1, r2, p), 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 (@(r2, r1) lder(r2, r1, p), 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), p);
    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), p);
    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

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
% 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.
%
% Revised 8/22/2018

function retval = lder(r2, r1, p)
    r3 = r1 + r2;
    r(1) = r1;
    r(2) = r2;
    r(3) = r3;

    for x = 1:3
        E1(x) = p.D1(x)*(exp(-2*p.B1(x)*(r(x) - p.r10(x))) - ...
                    2*exp(-p.B1(x)*(r(x) - p.r10(x))));
        E3(x,1) = p.D3(x)*(exp(-2*p.B3(x)*(r(x) - p.r30(x))) + ...
                    2*exp(-p.B3(x)*(r(x) - p.r30(x))));
        E3(x,2) = p.C(x)*(r(x) + p.A(x))*exp(-p.Si(x)*r(x));

        if (r(x) <= p.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*p.B1(i)*p.D1(i)*(exp(-2*p.B1(i)*(r(i) - ...
                        p.r10(i))) - exp(-p.B1(i)*(r(i) - p.r10(i))));

        if (r(i) <= p.R_1(i))
            dE3(i,i) = -2*p.B3(i)*p.D3(i)*(exp(-2*p.B3(i)*(r(i) - ...
                        p.r30(i))) + exp(-p.B3(i)*(r(i) - p.r30(i))));
        else
            dE3(i,i) = p.C(i)*exp(-p.Si(i)*r(i))*...
                        (1 - p.Si(i)*(r(i) + p.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
% 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.
%
% Revised 8/22/2018

function retval = rder(r1, r2, p)
    retval = 1 / lder(r2, r1, p);
end%function