Fukui matrix formation and diagonalization

Questions
#1

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

#2

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.

#3

Thanks for your response!

#4

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

#5

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.

#6

Thanks. I’ll give it a look.

#7

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

#8

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

#9

Thank you for the reference.

#10

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

#11

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)