Surface forces

[moritzgubler]

This is my implementation of the idea I had of calculating forces with surface integrals. It works quite well and the forces are more accurate. Also, mrchem struggles to compute forces when the world precision is 1e-7 or smaller. Then it requires a lot of memory (more than 32 GB for a H2 molecule). This is not the case with my approach . Geometry optimizations seem to be reliable if the stopping criterion is chosen 10 * world_prec.

I have already put all the parameters into the parser, so it is quite simple to test and run my code (look for the section "- name: Forces" in the template.yaml input parser file).

At the moment, there is two things that need to be improved:

1. I did not figure out a way to get access to the exchange correlation density and potential, so I had to change the xc_operator classes a bit in order to get access to them. It works now, but I did it a bit hacky.
2. I need the charge density, the electrostatic potential and the exchange correlation density and potential (the total one in case of an unrestricted calculation) to calculate the forces. A lot of the time when calculating the forces is spent recalculating these objects because I did not figure out a clean way to reuse them from the last SCF iteration. The time savings will probably not impact the entire SCF calculation significantly.
3. At the moment the method works only for LDA (this is also why some tests are failing). Extending the method to GGA would be super easy, the XC stress for GGAs is given here: Evaluation of exchange-correlation energy, potential, and stress Phys. Rev. B 64, 165110 – Published 4 October 2001 I am not sure how to calculate the exchange stress for hartree fock exchange. It shouldn't be too complicated but in a quick search I did not find a reference that derives it.

@stigrj Do you think my changes in the xc operator classes are fine.

[robertodr]

hei Moritz, could you add a quick description of how this differs from the Hellmann-Feynman calculation of the forces? Also, I have C++ code that:

• generates (at compile-time) the Lebedev quadrature points without reading them from a file,
• generates equal-weight grids on the sphere (Leopardi partitioning)

I can share both if it's of interest.

[moritzgubler]

I find i kind of usefule to be able to exchange the lebvedev file with something else and not have to recompile the code to test things. But if you like it more to avoid reading the file, I can also use your method.

What would be the purpose of the Leopardi points? are they more accurate than the lebvedev grids?

I sent you a small description of my method on zulip. I don't want to publish this yet.

[robertodr]

I find i kind of usefule to be able to exchange the lebvedev file with something else and not have to recompile the code to test things. But if you like it more to avoid reading the file, I can also use your method.

Yes, of course. The approach is the same, since the grids are known beforehand. If you tabulate all of them, there's not much you can swap them with, as their structure is not shared (that I know of) with other grids on the sphere. I dislike reading files at runtime. It's Lebedev, not lebvedev, by the way.

What would be the purpose of the Leopardi points? are they more accurate than the lebvedev grids?

In terms of accuracy they're similar. They're easier to generate though, so you're not limited to a pre-tabulated set of points and weights.

I sent you a small description of my method on zulip. I don't want to publish this yet.

Thanks!

[moritzgubler]

Can you send me the lebedev cpp library please you were talking about? I can also integrate that to avoid reading the file.

[robertodr]

here you go https://github.com/robertodr/balls

[moritzgubler]

Thanks, I integrated your library into my code.

[moritzgubler]

Open shell calculations now also work, I still need to implement the GGA stress density (see my question on zulip) and make sure, that the forces are calculated using the hellmann-feynman theorem in an exact exchange calculation.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 2.54735% with 1492 lines in your changes missing coverage. Please review.

Project coverage is 62.59%. Comparing base (df39ee2) to head (13b89d8). Report is 2 commits behind head on master.

Files	Patch %	Lines
src/surface_forces/SurfaceForce.cpp	0.00%	220 Missing :warning:
src/surface_forces/xcStress.cpp	0.00%	125 Missing :warning:
src/surface_forces/detail/lebedev_utils.cpp	0.00%	107 Missing :warning:
src/surface_forces/LebedevData.cpp	0.00%	75 Missing :warning:
src/surface_forces/lebedev.cpp	0.00%	36 Missing :warning:
src/surface_forces/detail/Lebedev_021_0170.hpp	0.00%	33 Missing :warning:
src/surface_forces/detail/Lebedev_023_0194.hpp	0.00%	33 Missing :warning:
src/surface_forces/detail/Lebedev_027_0266.hpp	0.00%	33 Missing :warning:
src/surface_forces/detail/Lebedev_035_0434.hpp	0.00%	33 Missing :warning:
src/surface_forces/detail/Lebedev_047_0770.hpp	0.00%	33 Missing :warning:
... and 31 more

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #488 +/- ## ========================================== - Coverage 68.66% 62.59% -6.08% ========================================== Files 194 233 +39 Lines 15287 16790 +1503 ========================================== + Hits 10497 10509 +12 - Misses 4790 6281 +1491 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[moritzgubler]

Hi @robertodr @ilfreddy I implemented now everything I wanted to, cleaned up to code and added a bit of documentation. The new method to calculate the forces now works for LDA and GGA functionals for both restricted and unrestricted calculations. The pulll request is now ready for review. Best, Moritz

[moritzgubler]

I previously tested GGA calculations only for H2 where the forces where the forces were correct. In more complicated systems this is not the case. (LDA is always correct). I will look into it and let you know once I have solved the problem.

[moritzgubler]

The GGA bug is fixed now, I forgot to add the divergence terms in the xc potential for the GGA case. They have been added now. MPI does not work anymore, it segfaults in the calculation of the kinetic stress. I suspect, that I try to access orbitals in an orbitalvector that are stored on another rank. @gitpeterwind Do you know what could be the problem?

[moritzgubler]

I can recreate the error with this input:

world_prec = 1.0e-4
world_unit = bohr

WaveFunction {
  method = blyp
  restricted = false
}

Properties {
  geometric_derivative = true
}
Forces {
    method = "surface_integrals"
    #method = 'hellmann_feynman'
    surface_integral_precision = medium
}

Molecule {
$coords
H 0.86387418 -0.02789966 0.13501906
Li -0.44463806 -0.04858572 -0.04959119
$end
}

I start the simulation with:

mrchem --dryrun hli.inp
mpirun -np 4 mrchem.x hli.json

[ilfreddy]

@moritzgubler would you mind adding also some testing to your code? According to the code coverage report, your patch is only marginally covered. Tests are essential to make sure the code does not get broken by others later on.

[moritzgubler]

I added a test, is there a way to look at the codecov report?