Implementation of DFT-D4 in the dftb+ infrastructure

[awvwgk]

• add DFT-D4 as dispersion correction to dftb+
  • add keywords in parser/input
  • Parameter setup for D4 calculator
  • coordination number and derivatives
  • Dispersion energy calculation and derivatives
    • two-body contribution
    • non-additive three-body contribution (ATM)
  • electronegativity equilibration charges (requires modified Ewald summation)
    • molecular
    • periodic
    • MPI parallel
• add detailed documentation on DFT-D4
• add tests for DFT-D4
• follow style guide

[aradi]

A few remarks at this early stage already:

• We'll have to see, whether D4 can be hidden beyond the usual dispersion interface. As it has a charge-dependence, we may need to handle it independently. (E.g. the Tkatchenko-type many body dispersion which is also being integrated has its own interface as well, as it is not compatible with the usual dipersion API).

• We use the NAG-compiler for testing purposes. However, that compiler unfortunately did not implement submodules yet, so for the moment, let's use conventional modules only.

[awvwgk]

All right, the submodule was anyway just temporary (to speed up compilation). I will remove it when I am done with the implementation.

D4(EEQ) can be hidden behind the usual dispersion interface, for GFN2-xTB the D4 needs a tighter integration to the SCC, but this is an issue for the GFN2-xTB implementation, not for DFT-D4.

[awvwgk]

@aradi I hope you are fine with compiling the DFT-D4 parameter file into dftb+. It not as heavy as the copyc6 from DFT-D3 but still 6k lines.

[aradi]

Do we have any other possibilities? 😉 I think we can live with it it, although one could think about a more dynamical solution, where the data is read from disk on demand. The only question were then, how to find the external data reliably during run (hard coded paths derived from the install prefix?)

[awvwgk]

I tested my implementation against the reference implementation and the results match so far. I still need to have a closer look on the strain derivatives and the Ewald summation to make sure everything works as expected.

I will start with adding tests and detailed documentation soon and will be happy to address style issues.

[aradi]

@awvwgk Sebastian, short question: are we supposed to look at it and consider it for merging or is this still considered to be work in progress?

[awvwgk]

@aradi, from the implementation side I am done with this feature. I still need to figure out how and where to add tests and docs, but was occupied otherwise this week, so I couldn't really look into it yet.

[awvwgk]

@aradi now also the docs and tests are done, at least I hope so. I am looking forward to hear your comments.

[aradi]

I am not through the entire PR yet (several parts of dispdftd4.F90 are still missing), but submit the partial review, so that we can start discussion.

Two additional general points:

• We should have a discussion, whether zeroing at allocation should be generally adapted or not. I can see advantages (stable behaviour on all systems, even if an algorithm is buggy) as well as disadvantages (debug features like initializing with NaN at allocation won't work any more, so finding bugs is more difficult).
• The manual should relate the input parameters with the symbols in the reference paper.

[awvwgk]

@aradi Thanks for the review.

I carried the habit of zeroing at allocation along from previous projects, since I usually zero variables anyway. I can remove all source statements from the code. In case it should be kept, I could use a fypp macro to switch between zero and NaN initialization.

[awvwgk]

I was waiting for the second half of the review, but didn't expect to get it patched in directly, so many thanks to you.

I will look through the code again next week and see what I can do to move this forward.

Should I add the D4 damping parameters already?

[aradi]

Sorry, I thought in case of stylistic changes it is easier to push directly as to describe them. I should have, however, probably made a PR against your branch instead of pushing directly, so that it is less invasive. Anyway, for me this is almost ready to go, apart of the S9 parameter, which does not pop up in the equations in the manual, but have to be set somehow in the input.

As for the D4 damping parameters, yes, adding them to the manual would make sense, for sure.

[awvwgk]

Not a problem at all, that's why I allowed “edits from maintainers” in the first place.

[aradi]

Oh, I see S9 is already in the manual. However, I think we should describe it more, as it is not clear from the first glance, that it should be either 0.0 or 1.0 depending whether 3 body terms should be included. Or do any intermediate values make sense? If not, maybe we should make it a logical (ThreeBodyTerms = Yes/No) and set internally 0.0. or 1.0.

[awvwgk]

Somebody might want to try 0.9 or 1.1 for the s9.
I already got a request to implement this option (in Orca, not DFTB+).

[awvwgk]

So there are two comments left. Than we should be ready to go.

1. Anything I should do about testing with MPI?
2. Should I remove calc%gc? Having the possibility to scale the hardnesses from the input is not a bad thing, but should be used very cautiously. I would prefer to leave it in since it is consistent with the dftd4 and xtb implementation.

[aradi]

As for the MPI, I updated the test diamond test to be big enough to be tested with MPI up to 4 procs. That works OK. Of course, the DFT-D4 part itself is not MPI-parallelized. At some point, we should think about it, as for a 64-diamond cell the D4-calculation takes order of magnitude longer than the rest...

Leave the charge steepness as freely adjustable parameter.

[awvwgk]

That's caused mainly by the ATM term, there are quite a lot triples. Usually it is better (meaning faster) to disable it with SQM methods.

[awvwgk]

@bhourahine, this is an artefact from the numerical differentiation. Just move the atom with the spurious gradient contribution a bit and check again, the artefact will be gone.

diff --git a/geo.gen.template b/geo.gen.template
index 79ea3e72..755002a0 100644
--- a/geo.gen.template
+++ b/geo.gen.template
@@ -2,7 +2,7 @@
  Ga As
  1 1 0.00 0.00 0.00
  2 2 0.17 0.25 0.25
- 3 1 0.50 0.62 0.00
+ 3 1 0.50 0.42 0.00
  4 2 0.75 0.75 0.32
  5 1 0.38 0.00 0.50
  6 2 0.75 0.10 0.75

[aradi]

@bhourahine could you try it with a stepsize of 1e-3 (instead of the 1e-5) with the critical geometry? I made the experience, that the agreement between numerical and analytical forces degrade for D4 if the stepsize goes below 1e-3 Ansgström. I am still not sure, though, where it comes from. My guess would be more on numerical instability / noise in the D4 calculation as problems with the finite difference calculation. (Without D4 I get agreement up to 1e-9 with a step size of 1e-5...).

[bhourahine]

@aradi I tried the larger step size as well, without a significant change in results.
@awvwgk I guess the question is on which side the difference comes from at the original geometry. We could try a more sophisticated numerical differentiator (higher order, extrapolation to limit, ...) is its not obvious.
One test I haven't done recently is to look what other dispersion model give for force errors.

[awvwgk]

@bhourahine With the modified geometry we arrive at:

...
Difference between obtained and reference forces:
       6.639819041809E-10      -6.138754409316E-10      -2.032428021020E-09
       2.279917616678E-09      -7.210614050290E-10       7.227968154555E-10
       1.907800563361E-09      -1.178456929263E-09       1.849387968288E-09
       2.067615105170E-09       1.170149276086E-09       8.489017752633E-10
      -1.817979954023E-09       2.551010652718E-10       9.998537414679E-10
       9.454846679885E-10      -5.153777543621E-10       6.474279896918E-10
      -2.305635508904E-10       1.453821686301E-09       1.277009735803E-09
      -7.402150663249E-10       3.377009366590E-10       1.139094374381E-09
Max diff in any force component:
       2.279917616678E-09
...
Difference between obtained and reference lattice derivatives:
       2.672237854558E-06      -1.809075546172E-07       4.497044225212E-07
      -8.412911196160E-08       1.705759536695E-06      -5.068280745984E-07
       3.350999484108E-07      -6.422317633419E-07       3.701807014343E-06
Max diff in any lattice derivatives component:
       3.701807014343E-06

We should be able to compare to the dftd4 standalone and the xtb implementation of D4, we can for molecular systems or completely neutral periodic systems (which did match perfectly), but we cannot for periodic systems with charges because the Coulomb matrix for the EEQ is not identical, which is due to the different setup of the Ewald summation (even with the same Ewald splitting parameter).

[awvwgk]

We should discuss this at the devel-meeting in detail. There is still the option to use the same Wigner–Seitz ansatz like in dftd4 or xtb, than we should be able to reproduce the exact same energy and gradient. But I'm not sure if this is really an issue, the geometry optimizations I tested so far did converge smoothly.

[awvwgk]

I cleaned up the branch and rebased it against master. Also, I updated the legacy make build system.

[bhourahine]

The fix for the old make system is around line 135 of prog/dftb+/make.common add in

DFTD4

L_INCS += -I$(ROOT)/external/dftd4refs

If everyone is happy, we can then merge.

[bhourahine]

prog/dftb+/make.common
even (markdown, humph)