Offset error in testing dataset

[dmilardovich]

Hi everyone,

I trained a GAP for SiO2. It seems to work reasonably well in MDs, but the errors I find in testing datasets are quite high. Moreover, there seems to be a fix offset in the energy estimations.

I was wondering if you could tell me what could be causing this. I already re-calculated the free energies and tried both using the e0 parameter and adding the single-atom configurations to the training dataset. The results are the same. Next are the results and the training command.

Thanks in advance!

os.system("gap_fit energy_parameter_name=energy force_parameter_name=forces at_file=training_dataset.xyz gap={distance_Nb order=2 cutoff=4.0 n_sparse=15 covariance_type=ard_se delta=4 theta_uniform=2.0 sparse_method=uniform compact_clusters=T : soap atom_sigma=0.5 l_max=4 n_max=8 cutoff=4.0 cutoff_transition_width=1.0 delta=0.4 covariance_type=dot_product n_sparse=1000 zeta=4} default_sigma={0.002 0.2 0.0 0.0} config_type_sigma={surface:0.002:0.02:0.0:0.0:cluster:0.01:0.1:0.0:0.0} e0={H:-6.4926:O:-83.9795:Si:-29.7169} sparse_jitter=1e-8 gp_file=GAP_48.xml")

[gabor1]

the picture didn't come through. as discussed before, the preferred way is to include the single atom in the training data set and NOT set an e0 on the command line, and also specify a config_type_sigma for the energy of the single atom as 0.0001 (or you can do it in that single atom frame of the xyz file by including energy_sigma=0.0001 on the second line). so pls resend energy force and virial parity plots.

[dmilardovich]

Hi Gabor, thank you for your answer and sorry for the very late reply. I re-trained including the single atoms in the training dataset (and NOT setting the e0 parameter).

I used the GAP to run MD and then re-calculated the energies and forces of the created atomic configurations with DFTB. The, I used DFTB to run the 'same' MD (i.e., same system and temperature, but initial velocities are of course different) and used the GAP to re-calculate the energies and forces of the atomic configurations. I am very surprised to see that the energy offset seems to depend on which potential (i.e., GAP or DFTB) runs the MD. Moreover, the offset is very small for the very first steps by it quickly grows and remains in what seems to be a fixed value for the rest of the MD run. Why is this happening? If this was a problem with the free energies, the offset would be fixed, right?

The next figure might be hard to read, but red = GAP energy and blue = DFTB energy. X-axis = step number and Y-axis = Energy [eV].

However, there is no offset in the forces. Now, if there is no offset in the force, I would think that the offset is only in the energy. If it is only in the energy, I would assume it is in the free energy. But, as we saw before, the offset cannot be produced by the free energies, as it is not fixed.

Moreover, I repeated this test for a different atomic system with more atoms and the offset changed (to a larger value).

Could you please give me a hint at what could be producing this energy offset?

This is the gap_fit command:

os.system("gap_fit energy_parameter_name=energy force_parameter_name=forces at_file=GAP_TRAIN_initial.xyz gap={distance_Nb order=2 cutoff=4.0 n_sparse=15 covariance_type=ard_se delta=4 theta_uniform=2.0 sparse_method=uniform compact_clusters=T : soap atom_sigma=0.5 l_max=4 n_max=8 cutoff=4.0 cutoff_transition_width=1.0 delta=0.4 covariance_type=dot_product n_sparse=1000 zeta=4} default_sigma={0.002 0.2 0.0 0.0} config_type_sigma={surface:0.002:0.02:0.0:0.0:cluster:0.01:0.1:0.0:0.0:single:0.0001:0.001:0.0:0.0} sparse_jitter=1e-8 gp_file=GAP_initial.xml")

Thanks a lot for your time! Best, Diego

[albapa]

What is free energy here? My question is whether the forces you are training with are derivatives of the energy values with respect to the atomic displacements. You can perhaps test this for the DFTB potential energy surface using finite displacements.

[gabor1]

Interesting Diego. Let's try to get to the bottom of it. Firstly, can you please verify that the prediction of the single atom energies is indeed around 0.1 meV for all three elements. Then next, can you please tighten the energy sigma to 0.001 to see what effect that has. what is your k-point convergence (are you using more than one k-point?) another interesting experiment is to REMOVE all but one energy value from the training set, this probes the possibility that your energies and forces are not quite consistent (what Albert is point out above). One possible reason for that is if you are running your electronic structure calculation at finite smearing, but giving energies extrapolated to zero smearing.

[dmilardovich]

Thank you Albert and Gabor for your answers. I will answer your questions next:

(1) Free energies.

Si (DFTB): -29.7169 [eV] Si (GAP): -29.7210 [eV] Error: 4.1 [meV]

O (DFTB): -83.9795 [eV] O (GAP): -83.97963 [eV] Error: 0.13 [meV]

H (DFTB): -6.4926 [eV] H (GAP): -6.49259 [eV] Error: 0.01 [meV]

Why could the error in the estimation of the free energies vary so much if they were all trained with the same energy sigma?

---

(2) DFTB forces vs. PES gradient: I run the test you proposed with an O-O dimer and I got these results:

Dimer A: O-O at 2.000 [A] --> DFTB Energy = -171.4913 [eV] Dimer B: O-O at 2.001 [A] --> DFTB Energy = -171,4841 [eV] Dimer C: O-O at 2.002 [A] --> DFTB Energy = -171.4768 [eV]

From dimer B and A, the PES gradient is 7.199 [eV/A] From dimer B and C, the PES gradient is 7.300 [eV/A]

So, when using a central difference approximation for the gradient (average of numerical left-sided AB and right-sided BC), the result is 7.250 [eV/A] --- The DFTB force in dimer B is 7.240 [eV/A].

---

(3) Extrapolation of energies: We use gamma-point calculations. There is no smearing. You can see next an example of the energy output of one of our DFTB calculations, where all energies are the same.

Total Energy: -285.4313919801 H -7766.9834 eV Extrapolated to 0K: -285.4313919801 H -7766.9834 eV Total Mermin free energy: -285.4313919801 H -7766.9834 eV Force related energy: -285.4313919801 H -7766.9834 eV

---

I assume from these results that: (A) The forces are consistent with the energies and (B) The error in the free energy of Si is too high and I should re-train the GAP with a lower energy sigma for the free-atoms. Please correct me if I'm wrong.

[albapa]

We use gamma-point calculations.

Depending on the database, this could be a problem. Energies, forces and virials of smaller cells will be less accurate. Not sure how much of a problem it is in practice - if you only use configurations of 100s of atoms, it may be not a big issue.

[gabor1]

what is your smallest cell?

[gabor1]

can you verify that al your single atom configs (it's weird to call them "free energy"!!!) have the correct config type set?

[gabor1]

your finite difference test shows a discrepancy of 0.01 eV/A. now this doesn't mean that this is the error in the analytical forces, but it does suggest a further test: retrain with all force sigmas at least 0.1 (your surface ones are smaller) and see if the energy prediction improves.

(this is after fixing the Si isolated atom energy prediction)

[gabor1]

another thing that's nice to plot (I am not yet sure it will help resolve the energy shift) is the series of dimer curves between all elements. we don't expect that to match the QM (in this case dftb) of course, but it's shape should be not so crazy. do it with absolute energies as well, rather than shifting them to zero.

[dmilardovich]

Hi, thanks again for your answers.

• Regarding the gamma-point calculations, the dataset is composed of atomic configurations with 224 - 272 atoms, so I assume this should not be a problem. However, could it be a problem for the dimers and single atom configs?

• Yes, all single atom configs have the right type set: "single".

• I will re-training with a lower energy sigma for the single atom configs. If that fixes the Si isolated atom energy prediction, I will re-train with all forces sigmas at least 0.1 eV/A. I will post the results when they are ready, including the dimer curves between all elements.

[gabor1]

Presumably you used a very large lattice (or non periodic boundary conditions) when you computed the single atom energies?

[gabor1]

Also, do a test where you remove all energy data other than the isolated atom energies. If this fixes the energy of isolated Si, then there is some data inconsistency in the energies (I don't yet know what could cause is)

[dmilardovich]

Hi, the lattice for the single atom configs is 20 A:

1 Lattice="20 0.000 0.000 0.000 20 0.000 0.000 0.000 20" Properties=species:S:1:pos:R:3:forces:R:3 config_type=single energy=-83.9795 pbc="T T T" O 10 10 10 0 0 0 1 Lattice="20 0.000 0.000 0.000 20 0.000 0.000 0.000 20" Properties=species:S:1:pos:R:3:forces:R:3 config_type=single energy=-29.7169 pbc="T T T" Si 10 10 10 0 0 0 1 Lattice="20 0.000 0.000 0.000 20 0.000 0.000 0.000 20" Properties=species:S:1:pos:R:3:forces:R:3 config_type=single energy=-6.4926 pbc="T T T" H 10 10 10 0 0 0

---

I have just repeated these calculations with different lattices:

Si (cell = 5 x 5 x 5 A) energy = -29.6988 [eV] Si (cell = 10 x 10 x 10 A) energy = -29.7169 [eV] Si (cell = 20 x 20 x 20 A) energy = -29.7169 [eV] Si (cell = 40 x 40 x 40 A) energy = -29.7169 [eV]

O (cell = 5 x 5 x 5 A) energy = -83.9767 [eV] O (cell = 10 x 10 x 10 A) energy = -83.9795 [eV] O (cell = 20 x 20 x 20 A) energy = -83.9795 [eV] O (cell = 40 x 40 x 40 A) energy = -83.9795 [eV]

H (cell = 5 x 5 x 5 A) energy = -6.4953 [eV] H (cell = 10 x 10 x 10 A) energy = -6.4926 [eV] H (cell = 20 x 20 x 20 A) energy = -6.4926 [eV] H (cell = 40 x 40 x 40 A) energy = -6.4926 [eV]

I used periodic boundary conditions.

Based on these results, we can assume that the single atom energy calculations are OK?

[gabor1]

yes, the single atom energies look fine

[dmilardovich]

Hi,

(1) I reduced the energy sigma for the single atoms configs to 0.00001 eV. The single atom energy errors are now: Si --> 0.03 meV. O --> 0.004 meV and H --> 0.0001 meV. The total energy offset was reduced but it is still very noticeable:

The error in the forces did not change:

---

(2) I increased the forces sigma to 0.1 eV/A for all configurations (keeping the energy sigma of the single atom configs to 0.00001 eV). The forces error increased slightly but it did not solve the energy offset. The dissociation energies look reasonable:

---

(3) I re-trained using only the energies and I still see the energy offset.

---

Could there be something wrong with the data? Is there something else I should check?

[gabor1]

The fact that your single atom sigma is now exceedingly tiny, and yet your Si prediction for the single atom is still 0.03 means that something else in your dataset is conflicting, it is preventing you from nailing the silicon single atom energy. your dimer curves are nice.

remind me - is your dataset homogeneous (apart from the single atom configs), i.e. does it come from the same kind of process? can you summarise the different kinds of config types that you have ? one could do "binary search", i.e. cut the data set in half (if it has multiple config types, and eliminate entire types), and find where the single atom energy snaps to sub-meV accuracy (forget about any other quality metric).

But before you do that, I am weirded out by your last plot the two disconnected force blobs in your energy-only fit. please identify what kind of force components those are, what element, and what configuration (x y or z component), etc etc. it may or may not be related to the other problems, but it's definitely abnormal.

[dmilardovich]

Hi, the training dataset is composed of single atoms (Si, O and H), dimers (Si-Si, Si-O and O-O) and Si surfaces in different levels of oxidation (i.e., including different number of O atoms - as shown in the figures in the first post).

The testing dataset is composed of Si surfaces with very few O atoms (as shown in the first figure of the first post).

---

I checked those force blobs and they all correspond to the Y component:

---

Moreover, they are the Y component of the H atoms:

These are the forces errors in the Y-direction for one random configuration of the testing dataset. However, the errors for all configurations look extremely similar. Black = H, yellow = Si and red = O.

---

These H atoms are fixed, they are simply at the bottom of the surfaces to pasivate them and keep them in place. I will remove these H atoms, re-calculate the energies and forces for the training and testing datasets and repeat the tests.

I am not sure if removing these H atoms will solve the problem. In any case, there is no need for them to be part of the training/testing datasets, as they are not part of the process we intend to simulate. They are simply there because these data was produced in a previous project (in which I guess passivating the bottom of the surfaces was required).

If that does not solve the problem, I will follow your advice and perform a binary search to find out in which data type is the problem. I will post the results when they are ready. Thanks again for your time, Gabor!

[gabor1]

Ok, very bizarre. Which way are your surfaces oriented?

[gabor1]

So one thing is that these H atoms are agonist all identical and so you have very little diversity in your H atom environments. That doesn't explain why only the Y component is broken

[dmilardovich]

Hi Gabor, please find next a figure to see the orientation of the surfaces (all of them are oriented in the same way). The axis are depicted in the bottom-left.

I believe there is a logical explanation for why only the Y component is broken:

• There is a single layer of H, which means that there is no H-H interaction in the Z-direction (only H-Si interaction).
• There is a very weak H-H interaction in the X-direction, as the distance between consecutive H atoms is 3.881 A.
• There is a much stronger H-H interaction in the Y-direction, as the distance between consecutive H atoms is 1.584 A.

---

I will create a new small dataset without the H atoms, re-train and re-test with it. I will post the results

[albapa]

Sorry if I am misinterpreting the situation. You say you are constraining the Hs - which force are you using in the training? If it includes the constraint, it will not be consistent with the total energy.

[dmilardovich]

Hi Albert,

The H atoms are constrained during the DFTB MD (by setting Driver --> VelocityVerlet --> MovedAtoms). After running the MDs, I subsampled the trajectory and re-calculated the energies and forces with DFTB again (I know this is an unnecessary step, but I did it to make sure that the DFTB setup is consistent for all my data and that the constraint applied to the H atoms did not affect the energies and forces). The forces I use in the training and those coming from the DFTB single point re-calculations, which means that this constraint in the MDs should not affect them. Is this correct?

Is there a better way to keep the Si surface in place without using H or any other extra element?

[Amitcuhp]

Hi,

(1) I reduced the energy sigma for the single atoms configs to 0.00001 eV. The single atom energy errors are now: Si --> 0.03 meV. O --> 0.004 meV and H --> 0.0001 meV. The total energy offset was reduced but it is still very noticeable:

The error in the forces did not change:

(2) I increased the forces sigma to 0.1 eV/A for all configurations (keeping the energy sigma of the single atom configs to 0.00001 eV). The forces error increased slightly but it did not solve the energy offset. The dissociation energies look reasonable:

(3) I re-trained using only the energies and I still see the energy offset.

Could there be something wrong with the data? Is there something else I should check?

Can you give me the sample file and jupiter file for understanding the things I am not able to learn all soap gap2b.?Please send it will be very helpful.

[gabor1]

Those tests are all fine - but we know where the problem is, it is to do with those H atoms whose forces are those separate blobs when you train with energies only.

[gabor1]

The energy only training gives you GAP forces that are approximate, but generally good in terms of sign and order of magnitude. So it's telling you that those Hs should have high y component forces, this is what's consistent with the energies in your training set. But the DFTB forces you are giving it for those Hs are small.

As another check, can you REMOVE those bottom Hs, recompute the DFTB forces, and retrain (both with energies only, and with energies and forces). Of course we'll figure out how to make it work with the Hs in place, I agree having them there is a good idea, but let's get a clear confirmation that that's where the (as yet unknown origin) problem is.