gbasis icon indicating copy to clipboard operation
gbasis copied to clipboard

Integrals Tutorial Notebook

Open maximilianvz opened this issue 1 year ago • 7 comments

Steps

  • [X] Write a good description of what the PR does.
  • [ ] ~Add tests for each unit of code added (e.g. function, class)~
  • [ ] ~Update documentation~
  • [X] Squash commits that can be grouped together
  • [ ] Rebase onto master

Description

This is a separated notebook corresponding to the Integrals section of the notebook from PR #142 (closed).

Type of Changes

Type
:bug: Bug fix
:sparkles: New feature
:hammer: Refactoring
:scroll: Docs

Related Issue

maximilianvz avatar Jan 10 '24 17:01 maximilianvz

@maximilianvz and @FarnazH I submitted the changes to the integrals related notebook. This tutorial is more abstract than others. I changed this to be the second tutorial. I am sorry that it took me longer than I thought. I will now be making the evaluation tutorial. It should be ready tomorrow.

Things to note: The calculated dipole moment example results are not the same as the reference (from the .fchk). I don't see where I made a mistake if you could, please, review it (if we are not able to fix it we may delete this case).

marco-2023 avatar Jan 14 '24 00:01 marco-2023

@marco-2023, the notebook looks pretty good to me. I pushed some formatting changes and things of that kind. I have three comments:

  1. In the text below, output and downloaded don't currently link to anything. Also, there appear to be some unfinished comments in the code block ("formchk file available from" and "output from Gaussian16 calculation available from") Do you intend to change this? image

  2. It would be best to ensure notation is consistent throughout the notebook. For example, in the screenshot below, you can see you make reference to $\chi_v$, but I don't believe this gets used anywhere subsequently. Also, it appears that we keep switching from using subscripts $\mu$ / $\nu$ and $a$ / $b$ for notating basis functions. An example of this can be seen in Section 1.2.2.1 versus Section 1.2.5. It would be nice if this was made more uniform. image

  3. I cannot figure out the problem with the dipole moment example, but I can keep trying to play around with it. @FarnazH, it would be great if you could take a look, too.

maximilianvz avatar Jan 15 '24 19:01 maximilianvz

@FarnazH @PaulWAyers @leila-pujal @maximilianvz Dear all: I went through this notebook one last time. Thanks to @leila-pujal I fixed some wrong prompts at the begining. I am still not able to get the proper results for the dipole moment example.

If I understood @PaulWAyers, the electric component of the dipole moment for the molecule would be: $$\mu_{x,y,z} = \sum_{i,j}\gamma_{ij}\int\phi_{i}(\tau)\hat{\mu}_{x,y,z}\phi_{j}(\tau)$$

I changed the example to use this formula and still I did not recovered the correct results. I tried comparing this with the example using grid and is this part that is not accurate. I don't know if there is a problem with the moment integrals (unlikely) or if I misunderstood something (very likely).

marco-2023 avatar Jan 20 '24 18:01 marco-2023

@marco-2023 can you check that the formula I gave you gives the same result in the AO and MO basis? It should not matter if you transform the AOs to the MO basis....

I will try to think about what is going wrong. We may have a bug....

(P.S. You are sure that we're using the same "center" for both the gbasis and grid moments?)

PaulWAyers avatar Jan 20 '24 21:01 PaulWAyers

hi @marco-2023 (and others). I used the moment integral recently to compute the dipole moment analytically so maybe I can help with this issue. I am looking at the example where you compute the dipole moment in your jupyter notebook. If the goal is to match the dipole moment printed in a fchk the formula I use is: $ \sum_{i}^{atoms} Z_i (x - X_C)^{c_x} (y - Y_C)^{c_y} (z - Z_C)^{c_z} - \sum_{\mu \nu} \gamma_{\mu \nu}\int \phi_\mu (\mathbf{r}) (x - X_C)^{c_x} (y - Y_C)^{c_y} (z - Z_C)^{c_z} \phi_\nu (\mathbf{r})$

I assume the nuclear dipole moment is the first term and the electric dipole moment is the second term. I attach below the code I use to reproduce the values printed in the fchk but I list the differences here. I think maybe here I have some advantage because I believe the molecule I am using its center of masses is placed at the origin. But other than that I think you are missing exponentiating by the orders of the momentum. I think it would be instead of mol_cm = mol.atmasses @ mol.atcoords / np.sum(mol.atmasses)

mol_cm_x = (mol.atmasses @ mol.atcoords / np.sum(mol.atmasses) )  ** np.array([1,0,0])[None, :]
mol_cm_y = (mol.atmasses @ mol.atcoords / np.sum(mol.atmasses) )  ** np.array([0,1,0])[None, :]
mol_cm_z = (mol.atmasses @ mol.atcoords / np.sum(mol.atmasses) )  ** np.array([1,0,1])[None, :]

Then other than that, I use as center [0,0,0] when I do the electronic term but that may not be correct for you. I have attached code and a fchk for you to try (change the txt extensions to py and fchk to use them). I am thinking this on the spot so maybe I said something wrong in my comment. If there is something that sounds wrong to you or my code or you have doubts let me know.

0009_CH3Clfieldx+064.txt dipole_moment_fchk.txt

leila-pujal avatar Jan 20 '24 21:01 leila-pujal

@PaulWAyers @leila-pujal @FarnazH Dear all, at last, I managed to find the source of the discrepancies in the example.

  1. Between the Grid moment integrals and Gbasis moment integrals, the source is the numerical error (if a very fine grid is used the difference dissapears).
  2. Between the Gaussian dipole moment and our dipole moment the difference is the reference point used. In Gaussian instead of the center of mass, the center of the atomic charges is employed. (I did not find any (Gaussian) documentation on this and thus I was lost for a while)

I pushed the correction. Everything else was right in the example. I am very sorry for the confusion and thank you for all of your feedback.

marco-2023 avatar Jan 21 '24 21:01 marco-2023

Good catch. I actually prefer the definition with the center-of-charge, which is what one would use in classical electromagnetism (in physics). However, when thinking about the interaction of a molecular multipole moment with an external field for the purpose rotational spectroscopy, the molecule rotates around its center of mass (not charge), ergo the "chemical" convention.

The good thing is that the code works :-)

PaulWAyers avatar Jan 21 '24 22:01 PaulWAyers

@maximilianvz and @marco-2023 , thanks again for working on making the tutorials. Similar to my message for the evals notebooks, I reviewed this a while ago and everything looks good to me. I made some minor changes to the location of the files and some formating. I will merge this once the checks pass. Thanks again.

leila-pujal avatar Apr 29 '24 15:04 leila-pujal