Fitting to ab initio data with forces


I’m fitting a potential for the cross LJ terms between the interaction of a gas molecules and a carbon structure.
I’ve obtained the interaction energies by performing DFT single point calculations moving the molecule towards the structure.

I successfully obtained the LJ parameters fitting only to energies and weighting the configurations with Boltzmann weights.
The interaction energies curves are pretty close to my DFT data, and I’m starting my fit with forces with these parameters.
However, when I perform the fit including the forces, the parameters go crazy (negative epsilons, sigmas going to zero etc…), even if I try to fit just epsilons, or just sigmas, or just a few parameters at a time.
I’m suspicious that I might be doing something wrong in the input of the atomic forces, or that the forces are not being weighted.

The way that I’m passing the forces to GULP is in the observables section, by adding a “gradients” section with the force on each atom.
I’m using the flags “fit conv” to perform the fit, as I don’t want the lattice parameters to be optimized.
I have something like:

fit kcal conv
50.00000000 50.00000000 22.10400000 90.00000000 90.00000000 90.00000000
fractional 1
C1 core 0.135599999 0.999299999 0.556731813 0.0000000 1.0000 0.0000
C1 core 0.132780000 0.027500000 0.556731813 0.0000000 1.0000 0.0000

C1 core 0.134779999 0.985120000 0.501176257 0.0000000 1.0000 0.0000
C2 core 0.000000000 0.000000000 0.000000000 0.0000000 1.0000 0.0000
O2 core 0.000000000 0.011499999 0.954985522 0.0000000 1.0000 0.0000
O2 core 0.000000000 0.988500000 0.045014477 0.0000000 1.0000 0.0000
gradients eV/Angs
1 0.133382 0.051353 -0.003565
2 0.105167 0.085753 -0.003581
3 0.085530 0.116774 -0.003956

-0.04897332 1.760000

then the file follows with several configurations, until

lenn epsilon zero 12 6
C1 core C2 core 0.41327703E-02 3.3085125 0.000 12.000 1 1
lenn epsilon zero 12 6
C1 core O2 core 0.10017456E-01 2.9127077 0.000 12.000 1 1

Any ideas on what I might be doing wrong?

Thanks in advance!

Hi Henrique,
There isn’t anything obviously wrong with your input & so it probably needs a careful analysis which is beyond the scope of a forum. Some general remarks:

  1. Make sure that your energies from DFT are binding energies with respect to the isolated atoms otherwise they won’t make sense for a short range potential.
  2. Fits leading to strange parameters usually mean either a) your observables are insufficient to constrain the fit; b) the observables are inconsistent with each other; or c) the potential model can’t fit all of the data and so it ends up in a bad place (sometimes you need to put more weight on the most important data near the minima).
    Hope that’s some help.

Hi Julian, thank you for your reply!

I’ve been experimenting with the fitting of this system, and found that when I weight the energies in the observables section (by adding real number after the energy) the forces does not seem to be weighted.
Is this indeed the case?

If it is the case, is there a way to weight the forces too?


Hi Julian,

Nevermind, I just found that I should add the force weight for each atom in the gradients section.
I’m going to try that, I’m just posting in case anyone else come across this problem.


Hi Henrique,
Weights can be set in a couple of ways:

  1. Change the default weight prior to the data being read in using the “default_weight” option.
  2. You can specify the weight for specific observables with the “weight” option.