MRChemSoft / mrchem

MultiResolution Chemistry
GNU Lesser General Public License v3.0
30 stars 21 forks source link

Surface forces (after rebasing #488) #497

Closed ilfreddy closed 2 days ago

ilfreddy commented 3 weeks ago

Taken from PR #488 from @moritzgubler and rebased with the AZORA functionality which has just been merged in master.

Rest of the description verbatim from #488

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).

It works with lda and gga functionals and for both closed and open shell systems. At the moment, there is one thing that might be improved in the future:

The exact exchange stress is not computed. Therefore forces can not be computed for exact exchange. It shouldn't be too complicated to implement, but in a quick search I did not find a reference that derives it.

codecov[bot] commented 2 weeks ago

Codecov Report

Attention: Patch coverage is 27.86561% with 1095 lines in your changes missing coverage. Please review.

Project coverage is 65.53%. Comparing base (764eb22) to head (dc430cf). Report is 1 commits behind head on master.

Files with missing lines Patch % Lines
src/surface_forces/xcStress.cpp 38.93% 80 Missing :warning:
src/surface_forces/detail/lebedev_utils.cpp 38.31% 66 Missing :warning:
src/surface_forces/LebedevData.cpp 13.33% 65 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_047_0770.hpp 0.00% 33 Missing :warning:
src/surface_forces/detail/Lebedev_059_1202.hpp 0.00% 33 Missing :warning:
src/surface_forces/detail/Lebedev_071_1730.hpp 0.00% 33 Missing :warning:
src/surface_forces/detail/Lebedev_083_2354.hpp 0.00% 33 Missing :warning:
... and 29 more
Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #497 +/- ## ========================================== - Coverage 69.16% 65.53% -3.64% ========================================== Files 202 241 +39 Lines 15646 17141 +1495 ========================================== + Hits 10822 11233 +411 - Misses 4824 5908 +1084 ```

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

moritzgubler commented 1 week ago

@ilfreddy I have found two minor bugs that I would like to fix:

  1. change the default value of the parameter radius_factor from 0.6 to 0.5.
  2. Force calculation fails when only a single atom is present in the system (The force is of course zero but it should nevertheless work). This is because the radius of the integration sphere (0.5 times the bond length) is not defined. The fix would be in the distanceToNearestNeighbour function which must return 1.0 or any other positive value when only a single particle is considered.

Should I create a pull request for this branch?

ilfreddy commented 1 week ago

@moritzgubler I think I have fixed what you asked. Could you check and the fixes, and then review the whole PR? Thereafter we can most likely merge it into master (assuming no test fails this time).

ilfreddy commented 1 week ago

@moritzgubler the test fails with the new default, so I explicitly used 0.6 in the test. Although the gradient should in principle be the same, I assume this is a low prec calculation with an SCF which is not converged either. I therefore assume it is not granted the computed gradient makes physical sense at this stage.

moritzgubler commented 1 week ago

In order to reduce testing time, I did not converge the SCF in the test (I basically copied the test from the Hellmann Feynman forces). The resulting forces do not make physical sense and even depend on the shape of the surface integral but I think this is ok for testing.

ilfreddy commented 1 week ago

If you are OK with my fixes and you can review the whole PR, I believe this gan go now, unless @stigrj wants to have a look first.

ilfreddy commented 2 days ago

If @moritzgubler and @stigrj are happy this one can finally go for my part.

moritzgubler commented 2 days ago

The changes look good to me.