mrchem
mrchem copied to clipboard
Surface forces
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:
- 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.
- 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.
- 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.