UT theoretical chemistry code forum

I'm working on a bit of a different situation from the solid state systems usually studied with the Bader analysis. We're simulating the solvated electron in liquid ammonia, where the PAW basis set of VASP is crucial to allowing the electron to fully delocalize. In addition, the spatial extent is very much of interest here, so I'd like to know more about the Bader volumes themselves, so I've been running the critical points calculation to try and get a better understanding of how it looks in space, but it's giving some errors.

The relevant excerpt from the output:

=========================================
Number of critical point count: 13958
Number of nuclear, bond, ring and cage critical point counts : 352 6107 5807 1692
ERROR: TWO TYPES OF CP ARE TOO CLOSE TO EACH OTHER.
The CPs 405 433 have number of negative eigenvalues: 0 1
STOPPING. NO MORE CP WILL BE OUTPUT.
After a round of reduction
Number of atoms: 258
Number of critical point count: 112
Number of nuclear, bond, ring and cage critical point counts : 66 17 23 6
WARNING: There is more nuclear CP than number of atoms!
ERROR: FAILED Poincare Hopf Rule and Morse Relationship
Writting critical point output files
Critical point information are written in file: CPFU00.dat
Network Graph is written in file: LINK00.dat
outputting debugging information to allcpPOSCAR
writting atomic locations
end of critpoint_find
Heuristic features not turned on. Exiting
=========================================

And for reference, the system has 258 atoms and 192 bonds. If nothing else, it's clearly failing to correctly account for the CP for some reason.

The thing catching my eye is that it's saying that "There is more nuclear CP than number of atoms", but that is not true? 66 nuclear critical points found and 258 atoms.

I'm also not sure why it's saying that it failed the Poincare Hopf Rule. My reading said that this theorem states that you can't create a closed volume around a saddle point without also including at least one minimum or maximum, and I'm assuming this is in reference to the Bader volumes being constructed? My system is weird, having a diffuse electron spread across most of the system, but I don't see why that would cause it to fail the Poincare Hopf Rule. When I looked at the code, it seems like this error gets thrown after simply having a large number (>100,000) for cptnum, which seems odd. I'm guessing that it's running the reduction of the CP and then checking if the number of CP is still high, which would indicate failure? But coupled with the code stating that there are more nuclear CP than atoms when that's not true, I wonder if some of the logic for this module is flawed. Or I could be wrong, but either way, could anyone help me get the critical point module working?

For what it's worth, the actual charges, the ACF and BCF, perform flawlessly. ACF has exactly 258 volumes, it computes reasonable charges for every atom. BCF has 260 volumes, which matches expectations for having 2 free electrons in the system. It only seems to be the critical point module that's having issues, but I'd really like to take a look at the critical points.

When we see the line "ERROR: TWO TYPES OF CP ARE TOO CLOSE TO EACH OTHER." it means two neighboring voxels both have a critical point and they are of different types (in this case, bond and nuclear).
So far I haven't seen a case when the code says this and the CHGCAR don't actually have two CPs adjacent to each other.
This is also very common, especially if the CHGCAR doesn't have a lot of grid points or if ENCUT isn't extremely high.

In some tricky systems, even ENCUT = 800 with NGXF set to 30 per angstrom isn't enough

Just to confirm my understanding. The code does it's run and finds, in this case, 13,958 critical points before finding 2 CP next to each other and stopping, not looking for anymore past that; and that's why when it runs the reduction afterwards, the final CP count is so low, and presumably why it's failing Poincare Hopf (like there "should" be another point there, but it stopped early and never found it, so it appears to be violating the rule)? As in, if it weren't for these 2 CP being next to each other (and lets assume they are the only problematic pair in the frame), then it should find the correct number of CP?

Does it just look for voxels directly adjacent? Scientific accuracy aside, could I "cheat" it by taking the CHGCAR and running a linear interpolation so that there's at least one voxel in between and then it would run?

And just theoretically, how is it that the charge analysis is running just fine when the CP is failing? If I don't have the resolution to resolve the CPs properly, shouldn't that imply that I don't have the resolution to properly run the charge analysis?

I think I can comment on the last point. The charge density analysis is pretty well tested in terms of finding Bader volumes and charges. And it is completely separate from the critical point code, which is done separately as a post-processing tool. The critical point code was written and tested by wc5879 (Ray) but has not be extensively tested by others. It is certainly possible that the critical point code can find spurious critical points while the standard Bader analysis is robust. In fact, one reason that our standard code is pretty robust is that we do not worry about specifically finding the boundary of Bader volumes or the critical points along them; instead, we just focus on which parts of space are inside each volume.

Anyway, we really appreciate your feedback and I hope that you can follow up with Ray to resolve this issue and maybe improve his code.

you can use "-ignore_cp_conflict" and "-autocp"

Like Graeme said, a CHGCAR not fine enough for CPs could be more than enough for accurate Bader charge and volume analysis.

The reason why Bader charge and volume are still accurate is because the CPs that cause troubles are typically in extremely low electron density regions (orders of magnitude smaller in charge density comparing to bond CP between atoms that corresponds to real chemical bonds, and some times only marginally higher than the default vacuum charge density).