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!
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
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.T nocc_alph = ccdata.homos + 1 mocoeffs_occ_alph = mocoeffs_alph[:, :nocc_alph] total_density = mocoeffs_occ_alph @ mocoeffs_occ_alph.T if unrestricted: mocoeffs_beta = ccdata.mocoeffs.T nocc_beta = ccdata.homos + 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)