Orbital coefficients

[lohy]

Hi,

Is it possible somehow to get the orbital coefficients out when calculating the eigenstates of a molecule?

I cannot find anywhere in the documentation that mentions this, so it might be that you call it something else.

thank you in advance,

Louise

[Anders Blom]

As you say, we might call it something else, but therefore we need to know how *you* define these coefficients :-)
And, it might help to know what you need them for, i.e. what analysis you can perform with them.

On the surface it sounds like you can just extract the Mulliken populations to know how much s, p, d contribution there is to a state, although it's a bit of a rough analysis.

[lohy]

Thank you for getting back to me. :-)
Say we have a C-O molecule. I define the orbital coefficient as how much the atomic orbitals of C is contributing to the molecular orbital(MO) - at e.g. the HOMO of the molecule. This I wanted to use to compare the HOMO of the molecule at different electric fields.

I also tried to use the 1D projector to compare the MOs(eigenstates), but got a bit suspicious about the unit (psi Å^-1.5).  Maybe you can clarify this?

I hope this makes sense.

[Anders Blom]

Lucky for you I actually wrote a script to do this about 10 years ago, but nobody seemed to need the functionality so we never made it a feature in the software. The usage is pretty obvious; the return is a list (one entry per eigenstate) of arrays, which contain the "contributions" (projections) onto the angular momenta s, p, d, and f. I tested with NH3 from the example in the manual (https://docs.quantumatk.com/manual/Types/ElectronDensity/ElectronDensity.html) and the run this slightly clumsy code

Code: python

from PartialNorm import getPartialNorms
molecule = nlread("nh3.hdf5", MoleculeConfiguration)[0]
# HOMO state is 3
ix = 3
partial_norms = getPartialNorms(molecule, projection_atoms=[0])
print(partial_norms[ix])
partial_norms = getPartialNorms(molecule, projection_atoms=[1])
print(partial_norms[ix])
partial_norms = getPartialNorms(molecule, projection_atoms=[2])
print(partial_norms[ix])
partial_norms = getPartialNorms(molecule, projection_atoms=[3])
print(partial_norms[ix])

We can then see that the HOMO level derives almost 90% from Nitrogen (mostly p-states), but if we change to ix=4 we note that LUMO has a lot more H contributions (and the LUMO state is, as expected, dominated by s-states). Hopefully this is helpful for the analysis of your molecule too. This function can be used for a lot of other cool things, like mapping out where s and p states dominate in the band structure of GaAs (RGB color scale mapped from spd contributions) as in the attached picture!

[lohy]

Wonderful thank you!  - you had a great idea 10 years ago 

Just to be sure, the script "only" needs the configuration and the energy state of the molecule I am interested in? and then I can choose which atoms I want to project on? where you have e.g. 0, can I then have a list or do I need to have a line for each atom in the molecule?

[Anders Blom]

It just need to have converged the molecular calculation. You can project on individual atoms (the example shows projection on N, H, H, H through the numbers 0,1,2,3.) or a collection of atoms by giving more than one index in the list.