thanks to the openKIM community I could finish my potential model.
I have verified that the forces, energies and particleEnergies are
correct by using LAMMPS. It produces exactly the same numbers as IMD,
so the values are correct.
But when I run the vc_forces_numer_deriv test provided at KCC I see
large deviations and the alpha at the end is of the order of 10^-3 and
not 10^-10 as expected. The forces that are computed differ considerably, e.g.:
Atom Type Dir Force_model Force_numer Force diff pred error weight stat
1 0 1 3.877238373442486E-01 -2.731087349161498E-01 6.60833E-01 5.55112E-17 1.22997E+15 FAIL
2 3.838638296887518E-01 -2.867688962737701E-01 6.70633E-01 1.66533E-16 1.29149E+15 FAIL
3 3.933035953364041E-01 -2.813365660120671E-01 6.74640E-01 2.03965E-12 1.37934E+11 FAIL
This is however not the case if I disable the embedding part of the potentials.
For pair potentials the test gives the expected results. When I enable the embedding part
it goes from alpha~10^-12 to 10^-3.
I guess that it might have something to do with the process_dEdr values.
Is there a way for me to check if I accounted for all contributions in this function?
The way I implemented this process_dEdr function is that after each force contribution I
also call the process_dEdr function with the value of the derivative. Do I have to take
care of something else?
This should not have anything to do with the process_dEdr function. The vc_forces_numer_deriv test only works with the `forces' array and `energy' value.
From the listing below I see that the numerically calcualted force is actually of the opposite sign from the analytically computed value. I wonder if there is a sign error on the embedding term. Did you try subtracting instead of adding the embedding contribution?
If the test does not depend on process_dEdr then I don't know what could be wrong.
As I said, the numerical values of LAMMPS and IMD are identical for forces, energy and particleEnergy.
I also checked this with an analytic potential by calculating these quantities by hand and they
are indeed correct.
Do you have any other ideas why this could be happening?
The numerical derivative code compares the forces computed by the model with a numerical derivative of the energy. It does not use the process_dEdr routines. The test is performed on a randomly perturbed FCC cluster where the atomic species are randomly selected from the species supported by the model. It could be that vc_forces_numer_deriv is probing your model in a region of configuration space (or with a species combinations) that is not visited by the IMD/LAMMPS test that you are doing.
Here is what I suggest you do:
1. The latest version of the vc_forces_numer_deriv code is attached. Please try that and see if you are still getting errors. (The new version of vc_forces_numer_deriv has some bug corrections and uses a better choice of the FCC lattice parameter size.)
2. If you are still getting errors, I suggest you modify the vc_forces_numer_deriv to write out the configuration that it is testing and run that through IMD or LAMMPS to see if you are getting the same forces. That should help pinpoint the problem.