Fitting EAM potential

Dear all.

I am trying to fit Nb EAM potential in GULP using this script:

fit conp opti prt_eam

cell
3.3008 3.3008 3.3008 90 90 90

fractional
Nb 0.0 0.0 0.0
Nb 0.5 0.5 0.5

observables
energy -15.14
elastic 3
1 1 246.6
1 2 133.2
4 4 28.1
end

buckingham
Nb Nb 1.0 0.1 1.0 0.0 4.5 1 1 1

manybody
Nb Nb 4.5

eam_functional square_root
Nb 1.0 1

# sum of cubic splines
eam_density cubic
Nb 0.1 3.0 1 1
Nb 0.1 4.0 1 1

Although the flags for parameters of the density and embedding functions are set, GULP does not optimize them (also prt_eam does not print any densities on atoms). Can someone, please, tell me what I am missing?

Thanks,

Roman

Hi Roman,
The issue here seems to be with the density and the use of the cubic form. If you use quadratic, as per the Finnis-Sinclair potential for Nb, then the energy is fine. On thinking about it, the issue is that for odd powers of (r-r0) then the potential will be negative and since you can’t have a negative density it gets put to zero. I’ll change the code for future versions to make the value the power of the absolute difference so that cubic is OK. For now this means that you’ll need to use quadratic or cubic forms (or place the cubic form on a grid). Sorry for that.
Best regards,
Julian