function f=fun(x,y,z,a,b,c)
% k=1e2;
k=20;
tmp = (c + cos(2*pi.*x) + a.*cos(2*pi.*y) + b.*cos(2*pi.*z));
tmp2 = exp(k.*tmp);
% dHda = -(k.*tmp2.*cos(2*pi.*y))./(tmp2+ 1).^2;
% dHdb = -(k.*tmp2.*cos(2*pi.*z))./(tmp2+ 1).^2;
dHdc = -(k.*tmp2)./(tmp2+ 1).^2;
f=dHdc;
end