# Fukui matrix formation and diagonalization

Is there a way to calculate and diagonalize the Fukui matrix in Q-Chem?

Q-Chem does not compute Fukui matrix (dP/dN, where P is the density matrix and N is the number of electrons) analytically, however it is possible to evaluate it numerically with a bit of post-processing as follows:
Run two subsequent jobs with a different number of electrons with `GUI=2` to save the density matrix in the formatted checkpoint file. Then use an external script to read in the density matrices and apply the finite difference formula (for example for a neutral system dP/dN ≈ 0.5*[P(+1) - P(-1)]) and diagonalize the resulting matrix. With Python this can be achieved using cclib and numpy.

Hi,

alternatively, you can run the jobs with `molden_format = true` for the different charges. Then you have to extract the section after `[Molden format]` to get the Molden files. If you have those you can use TheoDORE (https://theodore-qc.sourceforge.io/) to do the remaining parts using the `analyze_NOs.py` script.

-Felix

Note that the use of fraction electrons is also support (FRACTIONAL_ELECTRON in \$rem; see the manual). As popularized by Weitao Yang, these are used to construct E(N) plots to diagnose self-interaction error that manifests as deviations from linearity.

Thanks. I’ll give it a look.

Thank you for your reply. Is there a specific article you could suggest?

https://doi.org/10.1126/science.1158722

Thank you for the reference.

Does cclib natively know the format of the density matrix or do I need separate code to extract it?

Yes. Here is a starting point; I suggest passing the formatted checkpoint files created via `gui=2`, since the coefficients have higher precision:

``````from cclib.io import ccread

def calculate_total_density(ccdata):
unrestricted = len(ccdata.homos) == 2
# correct an annoying convention of parsing as [nmo, nbsf]...
mocoeffs_alph = ccdata.mocoeffs[0].T
nocc_alph = ccdata.homos[0] + 1
mocoeffs_occ_alph = mocoeffs_alph[:, :nocc_alph]
total_density = mocoeffs_occ_alph @ mocoeffs_occ_alph.T
if unrestricted:
mocoeffs_beta = ccdata.mocoeffs[1].T
nocc_beta = ccdata.homos[1] + 1
mocoeffs_occ_beta = mocoeffs_beta[:, :nocc_beta]
total_density += mocoeffs_occ_beta @ mocoeffs_occ_beta.T
else:
total_density *= 2

if __name__ == "__main__":
import argparse

parser = argparse.ArgumentParser()