help with writing a nebspline script in matlab

[egillsk]

I was trying to write a nebspline script in matlab (which is supposed to be a computer lab exercise in a chemistry class) but there seems to be something wrong. The spline is not continuos and takes very high and low values (on the order of 10 eV when the barrier is 0.5 eV and the forces are low). Something is wrong here. I both tried to do the cubic coefficients like in your papers and also as in the nebspline.pl script. However, I do not think those two formulas are the same though, and the results are quite different. Anyway, I have tried to make this work, but have not spot the error. Any help, asap, would be great so that I can finish this exercise for the students tomorrow...

here is the matlab code where I am trying to do it just like in the papers introducing this
G. Henkelman, G. Jóhannesson, and H. Jónsson, `Methods for Finding Saddle Points and Minimum Energy Paths', Progress on Theoretical Chemistry and Physics, ed. S. D. Schwartz (Kluwar Academic Publishers, 2000)
G. Henkelman and H. Jónsson, `Improved Tangent Estimate in the NEB Method for Finding Minimum Energy Paths and Saddle Points', J. Chem. Phys. Vol. 113, p. 9978 (2000):

clear all
close all

R = [0 0.9 1.9 2.5 3.8]; % the coordinates in Angstrom
E = [0 0.2 0.5 0.2 0]; % the energy of each image
%x = 0:0.1:10;
%f = spline(R,E,x);
F = [0 0.2 0.0 0.2 0]; % the parallel forces on each image
%F = zeros(size(R));

for i=1:length(R)-1
dR = R(i+1)-R(i);
F1 = F(i);
F2 = F(i+1);
U1 = E(i);
U2 = E(i+1);
d = U1;
c = -F1;
b = (3*U2 - U1)/dR^2 + (2*F1 + F2)/dR;
a = (2*U2 - U1)/dR^3 - (F1 + F2)/dR^2;
r = [R(i):0.1:R(i+1)];
e = a*r.^3 + b*r.^2 + c*r + d;
plot(r,e)
hold on
end



and here is the code when I am trying to mimic the perl script:

clear all
close all

R = [0 0.9 1.9 2.5 3.8]; % the coordinates in Angstrom
E = [0 0.2 0.5 0.2 0]; % the energy of each image
%x = 0:0.1:10;
%f = spline(R,E,x);
F = [0 0.2 0.0 0.2 0]; % the parallel forces on each image
%F = zeros(size(R));

for i=1:length(R)-1
dR = R(i+1)-R(i);
F1 = F(i)*dR;
F2 = F(i+1)*dR;
U1 = E(i);
U2 = E(i+1);
Fs = F1 + F2;
Ud = U2 - U1;
d = U1;
c = -F1;
b = 3*Ud + F1 + Fs;
a = -2*Ud - Fs;
r = [R(i):0.1:R(i+1)];
e = a*r.^3 + b*r.^2 + c*r + d;
plot(r,e,'r')
hold on
end


Best regards
egillsk

[graeme]

Perhaps I can help. There are some missing brackets in the equations in the NEB paper. The problem, and the correct equations are here: http://theory.cm.utexas.edu/henkelman/p ... _9978e.pdf

I have not run your script, but one thing that looks fishy to me is the forces. If this is a barrier, the projection of the force along the path should be first negative and then positive. Your data is all positive, which will lead to an unphysical wiggly path.

Hope this helps. If there are still problems, we can take a closer look. You can also create a small table with the distances, forces, and energy in a file called neb.dat and then call our spline script so that you can can compare directly to your implementation.

[egillsk]

Hi Graeme

Thanks a lot for your help. It finally worked. I also got some help from Andri through private emails. However, I think there is still a small error in the corrected formula in the Errata of the paper. But not in the perl script. If I changed the formula by having:
1. a negative sign in front of the 2 in the A constant
2. and by having (2F_{i}+F_{i+1})/R instead of 2(F_{i}+F_{i+1})/R, then it works.

Maybe you need to make a new Errata to the old Errata...hehe.. or tell me if I am wrong. Here is the new code below. I have tried to reproduce Fig. 5 in your J. Chem. Phys. (2000) paper. I also compare it with the cubic spline in matlab. The one with the DFT forces is similar as in your paper, but I guess my values which I read out from your figure are not completely right so this does not look like a converged NEB in my matlab figure. But I guess good enough.

____________________________________________
clear all
close all

R = [0 1 2 3 4 5 6 7 8]; % the coordinates in Angstrom
E = [0 0.1 0.4 1.8 3.3 3.4 2.4 1.99 1.96]; % the energy of each image
F = [0 -0.2 -1 -3.5 -2 1.3 1.8 0.2 0]; % the parallel forces on each image

for i=1:length(E)-1
for j = 0:0.01:(R(i+1)-R(i))
dR = R(i+1)-R(i);
F1 = F(i);
F2 = F(i+1);
U1 = E(i);
U2 = E(i+1);
d = U1;
c = -F1;
b = 3*(U2 - U1)/dR^2 + (2*F1 + F2)/dR;
a = -2*(U2 - U1)/dR^3 - (F1 + F2)/dR^2;
%r = [R(j):1:R(j+1)];
e = a*j.^3 + b*j.^2 + c*j + d;
coord = R(i) + j;
plot(coord,e,'r-')
hold on
end
end

z = 0:0.1:R(end);
hold on
y = spline(R,E,z);
plot(z,y,'k-',R,E,'o')
____________________________________________

best regards, egillsk

[graeme]

Great, thanks for the fix! You are probably the first and only person who has read that erratum (which was not actually published). I'll put up a fixed version.