DFT-D3 van der Waals does not work correctly with hydroxyl (-OH) groups

[mondracek]

The van der Waals forces in the FFvdW_?.xsf seem to be too small when calculated over an -OH group.

I first hit into this problem when trying to simulate a cyclodextrin molecule using the full-density based method. Trying to track the problem, I checked that even an isolated .OH radical (O+H pair of atoms) reproduces this problem. The resulting interaction of the CO tip with the hydroxyl group remains repulsive even for distances in which attractive interaction is expected, apparently because of a negligible (though non-zero) vdW component.

[NikoOinonen]

I assume this is using the DFT-D3 vdW, since you are talking about the full-density based model. I was also worrying about the magnitude of the vdW before in relation to the pyridine example. There specifically I found that the standard LJ vdW was too strong compared to the DFT-D3 vdW that was used in the original paper for the FDBM, which is what lead me to implement it also here. I think it might in general be the case that the DFT-D3 contribution is smaller than LJ.

If the vdW is too weak, there are several things you could try:

• It could be that instead of vdW being too weak, the repulsion is too strong. If you haven't tried yet, you could try adjusting the parameters of the Pauli integral, either increasing the exponent (small increments, like 0.05), or decreasing the prefactor.
• Try adjusting the parameters of the DFT-D3 (a1, a2, s6, s8): https://ppafm.readthedocs.io/en/latest/ppafm.ocl.html#ppafm.ocl.field.ForceField_LJC.add_dftd3
• Try switching back to the LJ vdW.

[mondracek]

@NikoOinonen thanks.

• As for your last suggestion, how do I do it? How can switch to LJ vdW with the full density in PPAFM? (Or perhaps it is actually DFT-D2 rather than LJ, as I definitely do not want the repulsive part of LJ to be combined with the Pauli repulsion from full density.)
• Regarding the suggestions to play with the parameters, well, no, I do not like it. I mean, I can try it just to find out more specifically where the problem is, but that's it. It is still a problem when the default set of parameters fails so badly. I can understand that of course I should not expect DFT-D3 to be absolutely perfect and fail-safe, like no method ever is, nor can any set of parameters be totally universal. One can always find some strange outlier molecule for which the method either fails or for which some tuning of parameters turns out to be the necessary and reasonable thing to do. However, hydroxyl groups are nothing that terribly unusual. If the default parameters fail for them, they are going to fail for most oxygen containing organic molecules. I really think there is a problem we need to solve.
• I will try to find out what the actual value of the C6 coefficient is for oxygen in the troublesome calculations and compare it to e.g. C6 for carbon (if there were order-of magnitude differences, that would indicate a problem). Besides, do you have any simulations with oxygen-containing organic molecule which you believe give reasonable results? I would like to check those if you do.

[NikoOinonen]

How can switch to LJ vdW with the full density in PPAFM?

Depends on whether you are using the CPU or the GPU version.

• CPU: Instead of ppafm-generate-dftd3, use ppafm-generate-ljff with the option --ffModel vdW and it will only output the vdW-part of the forcefield
• GPU: In AFMulator set fdbm_vdw_type = 'LJ'.

It is still a problem when the default set of parameters fails so badly.

Yes, probably at least the default Pauli parameters should be adjusted. But I'm not sure what they should be. In the papers where they used the FDBM, they mostly fitted the parameters to some DFT calculations, with quite some range in the optimal parameters. Currently, we have set them as prefactor=18, exponent=1.0, which is what the Perez group used as fixed values in their QUAM-AFM database.

[ondrejkrejci]

@NikoOinonen could you please specify DFT code, basis set and XC for the calculations for Fig. 2? @mondracek could you specify the same for your system? My personal guess is that the problem could be there.

[NikoOinonen]

@mondracek ran the DFT for the potential/density in Fig. 2, so I'm guessing they are the same.

I did not use the default Pauli parameters in that figure, however. There they were prefactor = 12, and exponent = 1.2.

[ondrejkrejci]

@mondracek ran the DFT for the potential/density in Fig. 2, so I'm guessing they are the same.

I did not use the default Pauli parameters in that figure, however. There they were prefactor = 12, and exponent = 1.2.

@NikoOinonen - So my guess is that you have used FHI-aims with PBE and Tkatchenko-Schefler for vdW. Isn't it? Would you be so kind to send the input files (control.in) to @mondracek , to see, if there are some differences in the options. Thank you!

[NikoOinonen]

@ondrejkrejci Yes, I used the PBE parameters for D3, but FHI-aims was not used in any part of the calculation. Martin ran all the systems with VASP I believe.

[mondracek]

• Guys (@ondrejkrejci and @NikoOinonen), I am confused now: Are we discussing Figure 2 of our paper (which seems to be alright) or the molecules for which I see some problems? Anyway, as for the DFT part, all my calculations have been done in VASP indeed. Which means plane wave basis set (cutoff energy 400 eV). The XC functional was PBE.
• The systems in which I observed problems were a cyclodextrin molecule at first and then, after I wanted to reduce the test to something very simple, a hydroxyl (.OH) radical. However, I have subsequently learned that there were actually two independent problems showing up, definitely for .OH and perhaps also for the cyclodextrin. One of the problems was the vdW-D3 interaction, the other was the DFT electron density from VASP, which sometimes ended up in a meta-stable solution, notably different from the sefl-consistent ground state. The wrong DFT density then translated into unrealistic Pauli repulsion. The latter problem is not a responsibility of our PP-AFM code, it entirely depends on the density that comes from VASP. I want to isolate the first issue only, that with vdW-D3. For that purpose, I have run simulations of a H2O molecule and I am going to upload the results, comparing forces obtained above the oxygen atom of a water molecule with three respective methods: LJ, FDBM + vdW-D2 and FDBM + vdW-D3.
• @NikoOinonen, I have now checked our paper and we claim the prefactor for the Pauli interaction to be A=1.2 there, instead of A=12 as you state above. Is it just a typo in the paper? If so, we should remember to correct it in the proof stage. By the way, I was using A=18 and exponent beta = 1.0 in my own simulations (of cyclodextrin etc.)

[mondracek]

Here is a figure with forces on H2O, comparing vdW-D2, vdW-D3 and Lennard-Jones. The full, dashed and dotted lines correspond to tip positions above the O atom of an H2O molecule with H atoms pointing down (in the xz-plane), above the O atom of a flat-lying (xy-plane) H2O molecule, and above an H atom of the same molecule, respectively.

I can upload the source data and input files for the simulations, if needed.

[NikoOinonen]

I have now checked our paper and we claim the prefactor for the Pauli interaction to be A=1.2 there, instead of A=12 as you state above

I see there is a typo in the figure caption. In the text it's correct, so A = 12.

[mondracek]

I am now trying to compare the forces fromppafm to DFT (enhanced with the standard vdW corrections as implemented in VASP). This issue is going to be more complicated then I imagined at the beginning. Anyway, I am not sure anymore there is literally a bug in our code. So I changed the label to possibly bug. What is clear so far is just that

• D2 and D3 vdW corrections are very different from each other
• D2 is reasonably similar to LJ
• D3 is very weak compared to D2 or LJ, but does not seem to be that different from Tkatchenko-Scheffler's correction as implemented in VASP. But this last point is currently a work in progress...

I need to run yet more DFT tests to understand what is going on. So, I am writing this comment to let you know what is the current status of the issue and that I am working on it, though there will be probably not much of contributions in the next few days, until I have the DFT results ready.

[yakutovicha]

@mondracek will check the issue again and close it if it is not relevant. Maybe add a small documentation section with images.

[mondracek]

Update to my previous comments:

• There is actually no vdW D2 correction implemented in PPAFM, only D3. So the D2 vs D3 comparison is moot. What I treated as D2 in PPAFM in my previous attempts was in fact the attractive part of the Lennard-Jones (LJ) interaction, but such a replacement of the D2 dispersion force model by LJ is not legitimate.
• D3 as implemented in PPAFM approximates both Tkatchenko-Scheffler correction and D2 correction implemented in VASP reasonably well. D3 is available only in the newer versions of VASP which are not available to me, so I cannot make a direct PPAFM vs VASP comparison for D3, unfortunately.
• There seems to be a systematic discrepancy between forces calculated with FDBM + vdW-D3 in PPAFM and DFT(PBE)-D3 in VASP in that the attractive forces are substantially weaker in PPAFM than they are in VASP This holds for several model systems in which the sample is represented by a free-standing atom or molecule. Nevertheless, the discrepancy does not originate from the van der Waals correction. If we wanted to address this, we should reopen the problem as a different issue.

Below, I will post a collection of figures showing the PPAFM vs VASP comparison of forces. Note: A=18 and B=1 were used as parameters of the Pauli repulsion in all the below graphs.

[mondracek]

Total forces including van der Waals for 6 different model systems, 4 molecules (benzene, bromomethane, carbon monoxide, and water) and 2 atoms (He, Xe). The tip is always a CO molecule, with the oxygen atom pointing towards the sample. On the molecules, the tip is placed over the C (in benzene), Br (in bromomethane), and O (in CO and H2O) atoms, respectively. The benzene and H2O molecules are flat-lying, CH3Br points with the Br atom towards the tip, CO is in an upright position wit the O atom pointing towards the tip.

[mondracek]

Van der Waals components only for comparison. This demonstrates that the D3 model of dispersion forces as implemented in PPAFM agrees reasonably well with the dispersion forces implemented in VASP.

[mondracek]

Forces without the van der Waals component. This shows the discrepancy between FDBM of PPAFM and VASP while also demonstrating that the problem is not with the dispersion (vdW) forces but somewhere else.

[mondracek]

This issue is ready to be closed. I will leave it hanging here open for a while so that there is a chance anyone notices the above figures. Question: Should I just leave the figures here or move them, for example, to a newly created Wiki page dedicated to them? I do not feel like these figures need to be part of the documentation but if you guys think otherwise, I can make that Wiki page. Comments?

[ondrejkrejci]

First of all - interesting that D32and TS do not match on CO-OC data - isn't there tilt of one of the CO molecules on D2?

What really surprises me is the huge difference between DFT and FDBM probably caused by charge (and multipole ) induction. I have no idea how to add these into the code easily. Probably some inspiration from classical force-fields or to try to go with some graph NN machine-learning force-field as an example in further future could be also an interesting project.

• I agree with closing this issue.
• I would really put these images on some dedicated wikipage as this could be really interesting.

[mondracek]

@ondrejkrejci Oops, there is a mistake in the CO-OC plot. I have put the van-der-Waals-only component there instead of the total force for DFT-D2. I will fix this. Edit: Done!

[ondrejkrejci]

Hi @mondracek - I have added all your images to forces description: https://github.com/Probe-Particle/ppafm/wiki/Forces#comparison-of-dft-vs-original-l-j-vs-fdbm-d3-implementation ; Thus we have all the important information both for us and also users, if they need it. Thank you for working on it!

[NikoOinonen]

I missed this since I was on holiday.

Good to see that the D3 is working more or less as expected, at least in the far range where it's the most relevant. And it actually should not completely agree in the close range, because in our implementation the probe particle does not count towards the bonding configuration of the sample, which it would in the DFT codes.

At the close repulsive range the FDBM seems to work pretty well, but I had no clue the discrepancy is so large in the middle range. In the original paper by Ellner et al. (https://pubs.acs.org/doi/pdf/10.1021/acsnano.8b08209) they show a good match, but conveniently it's only shown in the 2.8-3.4Å range. Still, I wonder if the match could be better with a different choice of the FDBM parameters.

[ondrejkrejci]

@NikoOinonen - yes, in here, we can see that already at 3.5 D2 and FDBM-D3 starts to overlay:

Both of those systems, are similar to what they had there. My personal guess is that they had a big discrepancies before and that is why they are not showing it, with motivation that bright lines are caused by repulsion and that the attractive forces are way weaker than the repulsive once. Thus we should be careful, where to fit those parameters ...

Onset of (even week) chemical interaction + inductional electrostatics (~1/r^4) are here to blame (my opinion).