Fix planewave based jellium model#1354
Conversation
The coulomb operator for jellium was hardcorded to the 3d case for planewaves: \sim 1/G^2, which gives incorrect results for 3d systems. The correct coulomb operator is now selected based on the dimensionality of the system.
The planewave basis set when used with an Ecut was truncated for the kinetic operator in the Hamiltonian, while the nuclear-electronic term was not truncated, and the electron-electron repulsion term was truncated incorrectly. Vee and Ven have been changed to match the kinetic planewave truncation.
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request introduces a generalized coulomb_potential_momentum function to compute the momentum-space Coulomb potential across 1D, 2D, and 3D systems, integrating it into various Jellium and plane-wave Hamiltonian models. It also adds support for filtering plane-wave operators within an energy cutoff. Feedback on the changes points out a performance bottleneck in the newly added filter_plane_wave_operator function, where repeated calculations and incremental additions to FermionOperator result in
The reduces the computational cost of filtering out high energy planewaves from an operator. Co-authored-by: gemini-code-assist[bot] <176961590+gemini-code-assist[bot]@users.noreply.github.com>
A nonexistent constructor was used when truncating the planewave basis of a FermionOperator.
mhucka
left a comment
mhucka
left a comment
There was a problem hiding this comment.
Thank you for doing this!
I have some comments about some possible slightly optimizations, but they might be better off left to another PR.
The only other thing is a request for formatting the formula in one of the docstrings. Other than that, it looks good!
| """ | ||
| if momenta_squared == 0: | ||
| return 0.0 | ||
| q = numpy.sqrt(momenta_squared) |
There was a problem hiding this comment.
Since momenta_squared is a Python float, if the value is a non-negative scalar, it should be possible to use Python's regular math.sqrt; this would be slightly faster as it will avoid some overhead in NumPy's version.
If this is changed to a math.sqrt, it also would be better to test that the value of momenta_squared is non-negative before taking the square root. The rest of this function tests other values and uses an else to deal with everything else, but it will be too late by then.
Finally, I see that the original code in this file used numpy functions and constants instead of Python math versions. I'm not sure if there was a reason for that – it might be a historical hold-over from pre-Python 3 days, or something else. If only the new code in this PR is changed to use the math versions, then that would be a bit inconsistent. So, an option might be to leave everything in this PR as you have it now (using numpy things), and leave the type changes to a different PR.
| return 0.0 | ||
| q = numpy.sqrt(momenta_squared) | ||
| if dimension == 3: | ||
| return 2.0 * numpy.pi / (volume * momenta_squared) |
There was a problem hiding this comment.
Similar as above, it should be possible to use Python's regular math.pi here; this would be slightly faster as it will avoid some overhead in NumPy's version.
| if dimension == 3: | ||
| return 2.0 * numpy.pi / (volume * momenta_squared) | ||
| elif dimension == 2: | ||
| return numpy.pi / (volume * q) |
There was a problem hiding this comment.
Ditto for the numpy.pi here.
| """Return the momentum space Coulomb potential for a given dimension. | ||
|
|
||
| For 1d systems, a soft coulomb potential is used with a regularization | ||
| parameter of a_1d: v(r) = 1/sqrt(r^2 + a_1d^2) |
There was a problem hiding this comment.
We should probably format this equation using LaTeX notation. This will render better when the docstrings are processed for the web documentation pages.
| ) -> float: | ||
| """Return the momentum space Coulomb potential for a given dimension. | ||
|
|
||
| For 1d systems, a soft coulomb potential is used with a regularization |
There was a problem hiding this comment.
Nit:
| For 1d systems, a soft coulomb potential is used with a regularization | |
| For 1-D systems, a soft coulomb potential is used with a regularization |
2 participants
This fixes #808, accounting for the 2 issues raised there: