will crash the Spyder and Jupyter kernels due to no_occupations. Is there another way to generate the natural orbitals for different molecules and find the interaction matrices (like through the .ao_overlap() in MintsHelper)? Thank you!
I can look into why you’re getting a “crash” with no further information, but two things I can say now:
I don’t know what “interaction matrices” are. Do you mean Hamiltonian matrix elements/integrals?
You can’t automatically determine the natural orbitals from setting up a molecular geometry and a basis set. You need to do some computation to determine the one-electron reduced density matrix and determine your natural orbitals from that. Different computations have different ways to get the density matrix back out. What method did you have in mind? (This is probably related to the crash - you haven’t computed NO occupations to return.)
Thank you for your reply! With regards to 1. I do mean the Hamiltonian matrix, and for 2. I’m not too sure how to approach getting the density matrix. I’ve used the .Da() function on wfn but haven’t been able to turn it to get it in array form. I understand that applying unitary operators to the density matrix to diagonalise it should produce the natural orbitals, is there a way to diagonalise the wfn.Da() density matrix for this purpose?
Your Hamiltonian integrals can indeed be gotten through similar functions as with ao_overlap. See the Psi4Numpy tutorials.
I don’t know what “haven’t been able to turn it to get it in array form” means, but I suspect you mean that it returns a Psi4 Matrix object. Again, see the Psi4Numpy tutorial I linked above, as well as the documentation on the NumPy interface.
As for the .Da() function, that returns the AO basis density matrix, but you want the MO basis density matrix. To get that, use .Da_subset("MO", "ALL").
I’ve tried using .Da_subset(“MO”, “ALL”) but get hit with the error
TypeError: Da_subset(): incompatible function arguments. The following argument types are supported:
1. (self: psi4.core.Wavefunction, arg0: str) -> psi4.core.Matrix
Invoked with: <psi4.core.Wavefunction object at 0x7fd940fc9890>, 'MO', 'ALL'
Which I assume means <psi4.core.Wavefunction object at 0x7fd940fc9890> is not the same class as psi4.core.Wavefunction, is there a way to change between the two if this the case?
No, the class is fine. The issue is the arguments. The supported argument types are a single string, not two strings, just as the error message says. This is my mistake. It’s just .Da_subset("MO"). There is no second argument to select only a subset of the MO density matrix.
In future, please enclose anything code-like in backticks (``) for formatting purposes. Single-backticks for in-line code, triple-backticks for a full block.
Ah I see, I’ve corrected that and seem to be getting density matrices that make sense. Would diagonalising this then produce the NOs? Also, is it also possible to produce and do this with a two-electron density matrix? Thank you so much again
Yes. The corresponding eigenvectors are the coefficients of the natural orbitals in the molecular orbital basis.
It isn’t clear to me what you mean by “do this” with a two-electron density matrix. If you want to obtain the 2RDM as a Python object, I would need to know which method you care about. It’s very uncommon to extract the 2RDM, and there’s no uniform interface for doing so. Even getting a uniform interface for 1RDM is a work in progress.
I’m trying to run some molecule simulations and the package in question requires the one and two-electron Hamiltonian matrices as inputs. I think in MintsHelper these are the sum of the ao_kinetic() and ao_potential() terms and the ao_eri() respectively and have used these for basis sets like 6-31G and cc-pvdz but was wondering if natural orbitals would be better. I’m not sure which methods exist for this but is it possible to get these Hamiltonian matrices from the 2RDM?
I’m still confused. I need more context about what your “big picture” problem is.
The molecular Hamiltonian is invariant to a unitary change of orbitals. The Hamiltonian using SCF orbitals is the same Hamiltonian as that using natural orbitals derived from some hermitian approximation to the 1RDM. If all you want is to feed the Hamiltonian to another package, you have no need to worry about the 1RDM or the 2RDM.
For approximate electronic structure methods that divide the orbitals into spaces (usually occupied and virtual), some choices of orbitals are better than others. If your package is making some assumption about “the first N orbitals are occupied”, then your electronic structure method isn’t just taking in a Hamiltonian. It’s also implicitly taking a partitioning of spaces. If this is what you’re talking about, tell me. Depending on what electronic structure method you’re worried about, you may also have to worry about orbital invariance.
While it is possible to get natural orbitals from the 2RDM, you do this by first determining the 1RDM from the 2RDM. This puts you in the same place as before.
Determining the Hamiltonian matrices themselves from the 2RDM of its ground-state wavefunction is an extremely difficult problem. Not all 2RDMs are even the ground-state wavefunction of any Hamiltonian, and finding the ones that are is a problem called pure v-representability, which has resisted most attempts to solve it. I’m not aware of any explicit solution for your problem (if this is your problem?), and the known explicit solutions are likely limited to toy systems.
Hi, thank you again. For the first point, is a unitary change of orbitals different to changing the basis set used? We have been inputting different one/two body matrices for different basis sets with corresponding changes in matrix size (for the one body Hamiltonians 6-31G’s size is 4x4, cc-pvdz is 10x10, and for pair natural orbitals is 2x2). Is changing between these different basis sets a unitary change?
For the second point, you’re right that the method does take an input N to choose the most occupied orbitals found (and will then use these N orbitals to construct the one and two body matrices). It then constructs qubit Hamiltonians based off of the HF orbitals and the N most occupied pair NOs.
For the last point, I’ve had a look at the source code and I don’t think the Hamiltonians are directly calculated from the 2RDM, but instead from the spin-orbital overlap integrals which are then used to generate the Hamiltonian. Are these overlap integrals more easily generated?
Changing between 6-31G and CC-PVDZ is not a unitary transformation. Those give different sets of atomic orbitals for a given molecular geometry, which are not related by any linear transformation, unitary or otherwise. You cannot write one set of atomic orbitals as a linear combination of the others.
It sounds like your method is truncating the space of virtual orbitals based on occupation numbers?
If this is a quantum computing project, then the Hamiltonians are certainly not directly calculated from the 2RDM. However, they cannot possibly be computed from the overlap matrix either. On its own, the overlap matrix gives zero information about the Hamiltonian. You need to combine that with OEI and TEI integrals in a not-necessarily orthogonal basis set for that to be at all useful. Look over your code again and figure out what it is you need. “How do I use the Psi interface,” I can help with. “What exactly do I need for my research project/external software,” I cannot.