Detect numerical failures due to the cavity early

[robertodr]

This PR provides a stopgap for cases where numerical artifacts in the cavity tesselation lead to problems, such as the one reported in #167. These are quite hard to predict given the molecular geometry, the radii set, and the average area of the finite elements. The strategy here adopted is to detect them and stop the calculation as soon as possible. This, as stated, is clearly a stopgap: the solution would be a more robust cavity generator.

Description

I have put two checks in place:

1. Once the S matrix is formed, we perform a robust Cholesky decomposition to check its positve-definiteness. We abort in case that's not the case.
2. Once the full PCM matrix (could be just S for CPCM) is formed, we compute the nuclear ASC and compare with the estimate given by Gauss' theorem. We abort for discrepancies in the absolute value of the difference larger than 10^-2.

How Has This Been Tested?

I have added two tests, to check that the failure conditions are correctly caught.

Questions

• Should the threshold for the Gauss' theorem test be more forgiving? Should I test on the absolute value of the relative difference instead?

Status

• [x] Ready to go
• [x] Cherry-pick to latest release branch

[MinazoBot]

	1 Warning
:warning:	Consider adding supporting documentation to this change. Documentation sources can be found in the doc directory.

Generated by :no_entry_sign: Danger

[robertodr]

Because I didn't refactor it... It can either go as a method of the Molecule class or as a free function in the utils folder. Any preference? I have no idea about the failures on Travis, but I'll investigate.

[ilfreddy]

No preference. I see Stig is lately moving a lot of functionality from classes to utils. So you can do that here too.

[robertodr]

Awesome! That was my preference too.

[robertodr]

What about my question on how to test the discrepancy? I've updated the CHANGELOG.md too

[codecov[bot]]

Codecov Report

Merging #170 into master will increase coverage by 2.51%. The diff coverage is 66.66%.

@@            Coverage Diff             @@
##           master     #170      +/-   ##
==========================================
+ Coverage   70.01%   72.52%   +2.51%     
==========================================
  Files          92       92              
  Lines        5592     5623      +31     
==========================================
+ Hits         3915     4078     +163     
+ Misses       1677     1545     -132

Impacted Files	Coverage Δ
src/utils/Molecule.hpp	100% <ø> (ø)	:arrow_up:
src/interface/Meddle.hpp	100% <ø> (ø)	:arrow_up:
src/bin/run_pcm.cpp	95.23% <100%> (ø)	:arrow_up:
src/bi_operators/IBoundaryIntegralOperator.cpp	100% <100%> (ø)	:arrow_up:
src/utils/Molecule.cpp	74.82% <100%> (+0.35%)	:arrow_up:
src/interface/Meddle.cpp	69.43% <45%> (-1.74%)	:arrow_down:
src/pedra/pedra_dlapack.f90	22.95% <0%> (+0.31%)	:arrow_up:
src/solver/CPCMSolver.cpp	74.19% <0%> (+6.45%)	:arrow_up:
src/utils/Atom.cpp	68.87% <0%> (+35.71%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update bb7e074...21ff187. Read the comment docs.

[ilfreddy]

I think we should test for the relative error (twice the charge, twice the charge error with the same cavity). The wonders of linearity ;-) An 10^-2 in relative error (1% missing charge) should be OK. Less than that is going to give many failures. And I am afraid this is still too tight

[robertodr]

@ilfreddy now implemented check based on relative difference. Please, have a final look-though