Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Docstring equations in evals and integrals modules #159

Merged
merged 15 commits into from
Jun 27, 2024

Conversation

gabrielasd
Copy link
Contributor

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)
  • Update documentation
  • Squash commits that can be grouped together
  • Rebase onto master

Type of Changes

Type
📜 Docs

Related

@gabrielasd
Copy link
Contributor Author

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.

@FarnazH FarnazH requested a review from maximilianvz January 12, 2024 23:18
@leila-pujal
Copy link
Collaborator

There are two TODO: equation in:

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.

@PaulWAyers
Copy link
Member

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.

@msricher
Copy link
Contributor

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.

@gabrielasd
Copy link
Contributor Author

gabrielasd commented Jan 15, 2024

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.

@msricher
Copy link
Contributor

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.

@gabrielasd
Copy link
Contributor Author

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

@maximilianvz
Copy link
Collaborator

maximilianvz commented Jan 15, 2024

@gabrielasd, all of your changes look good to me. I just have the following two requests before I approve it to be merged:

  1. This equation in contractions.py translates to:
Screenshot 2024-01-15 at 6 06 32 PM

      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.

  1. 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.

@gabrielasd
Copy link
Contributor Author

gabrielasd commented Jan 16, 2024

Hi @maximilianvz,
I've added the corrections. Please let me know if it is OK.

@maximilianvz
Copy link
Collaborator

Thanks, @gabrielasd. I'll approve this PR now.

maximilianvz
maximilianvz previously approved these changes Jan 16, 2024
@FarnazH FarnazH self-requested a review January 17, 2024 22:11
Copy link
Member

@FarnazH FarnazH left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gabrielasd thanks for your PR. Can you please address the small comments below?

Comment on lines 31 to 32
\phi_{\vec{a}, \vec{R}_A} (\mathbf{r}) =
N_{\phi} (\vec{a}, \vec{R}_A) \sum_i d_i g_i (\vec{r} | \vec{a}, \vec{R}_A)
\phi_{\vec{a}, \vec{R}_A, \mathbf{d}, \boldsymbol{\alpha}} (\mathbf{r}) =
N_{\phi} (\vec{a}, \vec{R}_A, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_i g_i (\vec{r} | \vec{a}, \vec{R}_A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. Both $\mathbf{r}$ and $\vec{a}$ are used in one euqation. Can you make everything $\mathbf{r}$, $\mathbf{R}$, and $\mathbf{a}$ in all equations?
  2. Instead of using subscripts for $\phi$, write it in the same style as $g_i$?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed here 7eb917b

Comment on lines 63 to 64
\phi_{j, \vec{a}, \vec{R}_A} (\mathbf{r}) =
N_{\phi} (\vec{a}, \vec{R}_A) \sum_i d_{ij} g_i (\vec{r} | \vec{a}, \vec{R}_A)
\phi_{j, \vec{a}, \vec{R}_A, \mathbf{d}, \boldsymbol{\alpha}} (\mathbf{r}) =
N_{\phi} (\vec{a}, \vec{R}_A, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_{ij} g_i (\vec{r} | \vec{a}, \vec{R}_A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please apply the same comments listed above.

Comment on lines 70 to 72
\rho(\mathbf{r}_n) = \sum_{ij} \gamma_{ij} \phi_i(\mathbf{r}_n) \phi_j(\mathbf{r}_n)

where :math:`\mathbf{r}_n` is the point at which the density is evaluated, :math:`\gamma_{ij}`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any reason to use $\mathbf{r}_n$ (with subscript n)? I recommend using just $\mathbf{r}$ to avoid confusion and be consistent. Please make the same change to all equations that follow.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think in the notes the notation was introduced in section 5.3 because of section 5.4, where it seems necessary for the notation of the derivatives of the density matrix (e.g. the equations for the stress tensor (Eq.33), the Ehrenfest force (Eq. 34) and others there).
However, I agree that for the docstrings in the code it may not be necessary, and definitely not for the density, and its 1st and 2nd derivative.
I'll change it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@PaulWAyers is the $r_n$ variable used in the equations for the positive kinetic energy density (Eq. 36), Ehrenfest force (Eq 34) Ehrenfest Hessian (Eq. 35) and stress tensor (Eq. 33), necessary? These are the equations also used in the docstrings of the corresponding functions in the code.
It seems what is being evaluated (an implemented) is the property for the trace of the density matrix $\gamma(\mathbf{r} , \mathbf{r})$, so it is not clear to me if the specification $r=r'=r_n$ is necessary.

Copy link
Contributor Author

@gabrielasd gabrielasd Jan 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made a partial fix making the change for the equations of the density and its derivatives, but I left out those for methods I asked Paul about (d345652)
For example for the positive definite KED:

t_+ (\mathbf{r}_n)

I was not sure if I had to turn the density matrix into the density function and the $\nabla_r \cdot \nabla_{r'}$ into $\nabla^2_r$ once I make this change

Comment on lines 205 to 229
.. math::
\begin{align}
&\frac{\partial^{L_x + L_y + L_z}}{\partial x^{L_x} \partial y^{L_y} \partial z^{L_z}}
\rho(\mathbf{r}_n)\\
&=
\sum_{l_x=0}^{L_x} \sum_{l_y=0}^{L_y} \sum_{l_z=0}^{L_z}
\binom{L_x}{l_x} \binom{L_y}{l_y} \binom{L_z}{l_z}
\sum_{ij} \gamma_{ij}
\frac{\partial^{l_x + l_y + l_z} \rho(\mathbf{r}_n)}{\partial x^{l_x} \partial y^{l_y} \partial z^{l_z}}
\frac{
\partial^{L_x + L_y + L_z - l_x - l_y - l_z} \rho(\mathbf{r}_n)
}{
\partial x^{L_x - l_x} \partial y^{L_y - l_y} \partial z^{L_z - l_z}
}
\end{align}

where :math:`L_x, L_y, L_z` are the orders of the derivative relative to the :math:`x, y, \text{and} z` components,
respectively.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. Change $\mathbf{r}_n$ to $\mathbf{r}$.
  2. The two lines seem to be numbered, like two equations.

Comment on lines 487 to 499
\begin{equation}
H[\rho(\mathbf{r}_n)]
=
\begin{bmatrix}
\frac{\partial^2}{\partial x^2} \rho(\mathbf{r}_n) &
\frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}_n) &
\frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n)\\\\
\frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}_n) &
\frac{\partial^2}{\partial y^2} \rho(\mathbf{r}_n)&
\frac{\partial^2}{\partial y \partial z} \rho(\mathbf{r}_n)\\\\
\frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n) &
\frac{\partial^2}{\partial z^2} \rho(\mathbf{r}_n)&
\frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n)\\
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variables in the derivatives are not correct. Please double check.

Copy link
Contributor Author

@gabrielasd gabrielasd Jan 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yup, sorry for missing this.
Addressed here: 8ca5e41

Comment on lines 154 to 155
\frac{\partial^{m_x + m_y + m_z}}{\partial x^{m_x} \partial y^{m_y} \partial z^{m_z}}
\phi(\mathbf{r}_n | \mathbf{R}_{A}, \mathbf{a}_j, \mathbf{d}_k, \boldsymbol{\alpha})
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can remove all other arguments of $\phi$ except $\mathbf{r}$ like when the basis functions appear in integrals and evaluations. It makes things simpler. What do you think?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made the change here 3fdabbc, together with the missing sum for the point charge.

@@ -271,6 +271,12 @@ def point_charge_integral(
):
r"""Return the point-charge interaction integrals of basis set in the given coordinate systems.

.. math::

\int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
Copy link
Member

@FarnazH FarnazH Jan 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gabrielasd, This is a matrix element for pairs of basis and the give point charge. The formula should be updated to better reflect that Returned value.

Copy link
Contributor Author

@gabrielasd gabrielasd Jan 30, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I understood the problem was a missing summation over the point charges.

I'm actually confused about the implementation of the formula for this function, point_charge_integral, Eq 47 page 21 in notes.pdf, in contrast with the one for nuclear_electron_attraction_integral, Eq. 47 page 21. To me the only difference in these equations is the value of the charge. But the implementations of both are different, with nuclear_electron_attraction_integral using point_charge_integral.

This function, point_charge_integral, like Farnaz points out, is returning a matrix of shape: number of basis functions by number of point charges, which is not reflected by the formula I added. On the other hand nuclear_electron_attraction_integral returns a square matrix of the dimensions of the basis functions.
point_chargeEq47
nuc_electron_attractionEq48

Copy link
Contributor Author

@gabrielasd gabrielasd Jan 31, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I ran the examples for both functions, and compared the shape of the returned variables; it seems to me that the documented return of point_charge_integral is missing one dimension:

eval_array : np.ndarray(K, N)

and should be np.ndarray(K, K, N) where K is the number of basis functions.

I'll add a comment to the docstring of Eq. 47 clarifying that what is returned are these (one-electron) integral matrix terms per every point charge, given that I'm not sure how to express it in the equation. For Eq. 48 I remember adding a sum over the charges in a previous commit.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @gabrielasd, I was checking your PR and specifically this part. I agree with you, I run the code and the point_charge integral function returns a matrix of elements being integrals of pairs of basis functions for each of the charges passed. Actually I was less sure about this line N is the number of coordinates at which the contractions are evaluated.` This would be charge points not coordinates right? Let me know what you think.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @leila-pujal, yes I agree with you.
It does seem a bit ambiguous since in this case the number of coordinates should match the number of charge points.

FarnazH added a commit to maximilianvz/gbasis that referenced this pull request Jan 18, 2024
Issue theochem#159 includes the formulas in the docstring, and the notebooks contain link to
the docstring, so there is no need to repeat these.
leila-pujal pushed a commit to maximilianvz/gbasis that referenced this pull request Apr 29, 2024
Issue theochem#159 includes the formulas in the docstring, and the notebooks contain link to
the docstring, so there is no need to repeat these.
Copy link
Collaborator

@leila-pujal leila-pujal left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I went over the changes added to the docstring in this PR. Looks like what was left to clarify was the $r=r'=r_n$ and the point_charge/nuclear_electron_attraction_integral documentation. I left a comment for the latter. For the $r=r'=r_n$ 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. The last thing I wanted to mention is that in some docstrings you start the math with \begin . I think that is only needed when you are writing on Latex and you need to activate the math mode to keep track of the equation but don't think we need this here. If you address these issues I think I will be able to merge the PR.

@@ -271,6 +271,12 @@ def point_charge_integral(
):
r"""Return the point-charge interaction integrals of basis set in the given coordinate systems.

.. math::

\int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @gabrielasd, I was checking your PR and specifically this part. I agree with you, I run the code and the point_charge integral function returns a matrix of elements being integrals of pairs of basis functions for each of the charges passed. Actually I was less sure about this line N is the number of coordinates at which the contractions are evaluated.` This would be charge points not coordinates right? Let me know what you think.

@gabrielasd
Copy link
Contributor Author

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?

@leila-pujal
Copy link
Collaborator

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?

@gabrielasd
Copy link
Contributor Author

Hi @leila-pujal, yes I think that is a good solution. I'll update the equations and get back to you.

@gabrielasd
Copy link
Contributor Author

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.

Copy link
Member

@PaulWAyers PaulWAyers left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The formula in @gabrielasd 's comment would not be right, so I don't think we want that formula. The current formula in line 624 looks right to me.

@leila-pujal
Copy link
Collaborator

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.

Copy link
Member

@PaulWAyers PaulWAyers left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The positive-definite kinetic energy corresponds to (one half of)
orders_one = [1,0,0] orders_two = [1,0,0]

  • orders_one = [0,1,0] orders_two = [0,1,0]
  • orders_one = [0,0,1] orders_two = [0,0,2]

This is equal, mathematically, to

$$ \iint \delta(\mathbf{r}_1 - \mathbf{r}_2) \left[ \nabla_1 \cdot \nabla_2 \gamma(\mathbf{r}_1, \mathbf{r}_2) \right] d \mathbf{r}_1 d\mathbf{r}_2 $$

where the Dirac delta function in the integral is usually replaced by just denoting that after evaluating the integral one sets $\mathbf{r}_1 = \mathbf{r}_2$.

I may have read it wrong but that seems to be the positive-definite kinetic energy implemented. If it were not, I would have expected that a test that the pos. def. k.e. density is nonnegative and integrates to the global k.e. would fail.

Changes:
- Add coeffs and exponent variables to the docstring equations defining
the generalized contraction functions in GeneralizedContractionShell based
on Eq. (11) from the PDF notes.
- Include the equation for the arbitrary order derivative of a basis function in
the docs of `evaluate_deriv_basis` (Equation 22 from the notes).
Changes:
- Add sum over point charges in point_charge_integral.
- Remove extra variables in evaluate_deriv_basis.
Notes:
- Change not applied to the docs of methods implementing derivatives of the
density matrix like posdef_kinetic_energy_density or those in stress_tensor.py
@leila-pujal leila-pujal force-pushed the docstring_equations branch from 9ea99ce to 3a885f2 Compare June 27, 2024 15:58
@leila-pujal
Copy link
Collaborator

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.

@leila-pujal leila-pujal merged commit 77bbbd4 into theochem:master Jun 27, 2024
0 of 10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants