theochem
Docstring equations in evals and integrals modules
When missing add equations to the docstrings of the high-level functions in the modules evals and integrals. The formulas were taken from the PDF documentation notes.
Checklist
- [ ] Write a good description of what the PR does.
- [ ] Add tests for each unit of code added (e.g. function, class)
- [x] Update documentation
- [ ] Squash commits that can be grouped together
- [ ] Rebase onto master
Type of Changes
| Type | |
|---|---|
| ✓ | :scroll: Docs |
Related
Hi @msricher and @leila-pujal although I added the equations to the docstrings from the notes.pdf I had doubts about the correspondence between the formulas and the examples (or the function arguments) in these two cases:
- nuclear_electron_attraction_integral (Section 6.4.1): I read the integral as referring to the interaction with one nuclei, but actual evaluation function, and the example indicates it is applicable to an array of nuclear charges and coordinates. Should I add a sum over nuclei?
- point_charge_integral (section 6.4): Should there be a charge variable there? The formula seems to me like the point charge is unitary, but the example mentions specific charge values.
The last doubt was about a minus sign in the electrostatic potential formula, Eq. 49, section 6.4.2. What is the side in the equation that is correct? On the right-hand-side of the equation the first and second term have opposite sign, but in the left-hand-side, the same two terms between parenthesis have the same sign.
There are two TODO: equation in:
- integrals/_diff_operator_int.py
_compute_differential_operator_integrals_intermediate - integrals/_moment_int.py
_compute_multipole_moment_integrals_intermediate
I have been looking more in detail for the angular momentum and the first one I think _compute_differential_operator_integrals_intermediate refers to $ (d_x)^{k+1}_{a_x b_x} = 2 \alpha d^k_{a_x+1, b_x} - a_x d^k_{a_x-1, b_x}$ (eq.52 Paul's notes). For _compute_multipole_moment_integrals_intermediate it looks it most likely corresponds to a specific $s_{ij}^e$ but I haven't looked at this one in detail. @PaulWAyers @FarnazH, do you think these equations need to be added to the docstring? @gabrielasd already took care of the main equations (e.g momentum integral, angular momentum integral etc) but these are intermediates we use to compute these integrals through the recursions. Since we are not including the recursions in the paper I wanted to check what you think about this.
More documentation, and more specific documentation, is always better but we shouldn't get too bogged down IMO since it's impractical to ever perfectly document anything.
Hi @msricher and @leila-pujal although I added the equations to the docstrings from the notes.pdf I had doubts about the correspondence between the formulas and the examples (or the function arguments) in these two cases:
- nuclear_electron_attraction_integral (Section 6.4.1): I read the integral as referring to the interaction with one nuclei, but actual evaluation function, and the example indicates it is applicable to an array of nuclear charges and coordinates. Should I add a sum over nuclei?
Yes it's a sum over nuclei.
- point_charge_integral (section 6.4): Should there be a charge variable there? The formula seems to me like the point charge is unitary, but the example mentions specific charge values.
The actual charge can be moved outside the equation, so it's trivial. In Libcint, for example, I implement it using just the integral of 1/|r - r_c|, (instead of -q/|r - r_c|), and then I multiply the entire result by -q. There is a negative sign because the interaction of an electron (-1 charge) with a positive point charge is attractive.
The last doubt was about a minus sign in the electrostatic potential formula, Eq. 49, section 6.4.2. What is the side in the equation that is correct? On the right-hand-side of the equation the first and second term have opposite sign, but in the left-hand-side, the same two terms between parenthesis have the same sign.
I think both sides of that equation are correct? If you bring all the negative signs to the outside of each term, then you have +Term1 - Term2 = +Term1 - Term2.
Thanks @msricher
- point_charge_integral (section 6.4): Should there be a charge variable there? The formula seems to me like the point charge is unitary, but the example mentions specific charge values.
The actual charge can be moved outside the equation, so it's trivial. In Libcint, for example, I implement it using just the integral of 1/|r - r_c|, (instead of -q/|r - r_c|), and then I multiply the entire result by -q. There is a negative sign because the interaction of an electron (-1 charge) with a positive point charge is attractive.
Got it!
The last doubt was about a minus sign in the electrostatic potential formula, Eq. 49, section 6.4.2. What is the side in the equation that is correct? On the right-hand-side of the equation the first and second term have opposite sign, but in the left-hand-side, the same two terms between parenthesis have the same sign.
I think both sides of that equation are correct? If you bring all the negative signs to the outside of each term, then you have +Term1 - Term2 = +Term1 - Term2.
Then shouldn't Term2 on the right-hand-side have a -1 on the numerator? like it had on the left-hand-side.
Then shouldn't Term2 on the right-hand-side have a -1 on the numerator? like it had on the left-hand-side.
Oh, you're right! I'd think both terms should end up being positive. Page 362 of Helgaker (eq. 9.7.4) gives ESP only in terms of the density $\rho(r_p)$, but I think this also implies that the second term should be positive. I'm not 100% sure though.
Then shouldn't Term2 on the right-hand-side have a -1 on the numerator? like it had on the left-hand-side.
Oh, you're right! I'd think both terms should end up being positive. Page 362 of Helgaker (eq. 9.7.4) gives ESP only in terms of the density ρ(rp), but I think this also implies that the second term should be positive. I'm not 100% sure though.
Yup, you are right. And based on the implementation, it is like that (I should have checked this before); Term2 is the variable hartree_potential on L128 and the function returns -(Tern1 + Term2).
So, I'll check that in the docsting I have the right signs for the right-hand-side of Equation 49 in the notes.
Thanks a lot @msricher
@gabrielasd, all of your changes look good to me. I just have the following two requests before I approve it to be merged:
- This equation in
contractions.pytranslates to:
This looks similar to Equation 11 from notes.pdf, but seems to be missing some variables (e.g., contraction coefficients, $d$, and the $\alpha$). This isn't your fault, because you haven't touched this file, but I would appreciate it if you could adjust the docstring so that it matches Equation 11 in the notes.
- I think it would be good to have Equation 22 from the notes in the docstring for
eval_deriv_basis, because this is the function that is called by all the higher level functions that require derivatives.
Hi @maximilianvz, I've added the corrections. Please let me know if it is OK.
For the r=r′=rn and your question about if by only using r in the notation, I think for the higher API functions we should change things like what you mentioned here. Right now I was thinking and I think the user by using the lower level functions can still evaluate the matrix and not only the trace.
Hi @leila-pujal, sorry, I still don't understand how to go about this fix. Do you mean that what I should do is turn the density matrix into the density function and the $\nabla_r \cdot \nabla_{r'}$ into $\nabla^2_r$ for the higher API?
Hi @gabrielasd, yes I meant that. Now that I read my comment again it is not entirely clear that I was pointing at this sorry. What do you think about this fix?
Hi @leila-pujal, yes I think that is a good solution. I'll update the equations and get back to you.
Hi @leila-pujal, I think I addressed the requested changes, but please let me know in case I missed something.
For the formula of the positive definite KED, in the end I didn't modified it. Based on the way it was implemented, I found that the equation in the notes made sense rather than making the change to: $t_+ (\mathbf{r})= \frac{1}{2} \nabla^2 \gamma(\mathbf{r}, \mathbf{r})$ but I can change to this one if it is preferred.
Hi @gabrielasd and @PaulWAyers thanks for catching this and sorry for the confusion. Just to double-check, equation in @gabrielasd message is wrong because that would refer to the laplacian of the electron density? The way it is implemented is that we compute for basis functions a and b derivatives (a)100/(b)100, (a)010/(b)010, (a) 001/(b)001 which is not the laplacian but the gradient of orbital a multiplied by the gradient of orbital b. Maybe I understood it wrong. Please let me know if this is not the reason why the equation in the message is wrong.
Thanks a lot @PaulWAyers for your answer and sorry for not getting back to you until now. I have forced push my branch that I rebase locally so we can merge it. Thank you so much @gabrielasd for working on this.