Surface forces (after rebasing #488) #497
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ 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 ☔ View full report in Codecov by Sentry. |
|
@ilfreddy I have found two minor bugs that I would like to fix:
Should I create a pull request for this branch? |
|
@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). |
|
@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. |
|
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. |
|
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. |
doc/users/user_ref.rst
Outdated
| :radius_factor: Sets the radius of the surface used in the computation of forces. The radius is given by this factor times the distance to the neariest neighbour. Must be between 0.1 and 0.9. This should rarely need to be changed. Different values can change the accuracy of the forces. | ||
|
|
||
| **Type** ``float`` | ||
|
|
||
| **Default** ``0.6`` |
There was a problem hiding this comment.
Didn't you change that to 0.5?
src/CMakeLists.txt
Outdated
| LEBVEDEV_SOURCE_DIR="${LEBVEDEV_SOURCE_DIR}" | ||
| LEBVEDEV_INSTALL_DIR="${LEBVEDEV_INSTALL_DIR}" |
There was a problem hiding this comment.
I think you can delete these lines
src/driver.cpp
Outdated
| @@ -57,6 +57,7 @@ | |||
| #include "qmoperators/one_electron/NuclearOperator.h" | |||
| #include "qmoperators/one_electron/ZoraOperator.h" | |||
| #include "qmoperators/one_electron/AZoraPotential.h" | |||
| #include "qmoperators/one_electron/NablaOperator.h" | |||
src/driver.cpp
Outdated
| #include <fstream> | ||
| #include <filesystem> |
|
|
||
| Forces { | ||
| method = 'surface_integrals' | ||
| radius_factor = 0.6 |
There was a problem hiding this comment.
when the radius_factor is set here, the test should still work when the default is changed to 0.5
There was a problem hiding this comment.
Indeed it does. With the added bonus that we test setting a non-default value 😄
Default value updated in the manual
python/template.yml
Outdated
| - name: radius_factor | ||
| type: float | ||
| default: 0.6 | ||
| docstring: | | ||
| Sets the radius of the surface used in the computation of forces. |
There was a problem hiding this comment.
Here, the default is also not 0.5 did you run the parameter update script?
There was a problem hiding this comment.
I fixed that too now :-) Thanks for pointing that out!
|
If @moritzgubler and @stigrj are happy this one can finally go for my part. |
|
The changes look good to me. |
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.