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.

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 analyze_NOs.py script.


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
        total_density *= 2
    return total_density

if __name__ == "__main__":
    import argparse

    parser = argparse.ArgumentParser()
    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)