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.

Thanks for your response!

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?

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
return total_density
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("outputfilename_ref")
parser.add_argument("outputfilename_plus")
parser.add_argument("outputfilename_minus")
args = parser.parse_args()
data_ref = ccread(args.outputfilename_ref)
data_plus = ccread(args.outputfilename_plus)
data_minus = ccread(args.outputfilename_minus)
density_ref = calculate_total_density(data_ref)
density_plus = calculate_total_density(data_plus)
density_minus = calculate_total_density(data_minus)
```