EzReson: An efficient program for chemical resonance analysis
The program provides a general and efficient tool to perform a resonance analysis of molecules. It translates the wave function obtained from a standard DFT or HartreeFock calculation to the chemically more intepretable Lewis structures.
Download link for the latest v2.0.1 version: https://github.com/yangwangmadrid/EzReson/releases/download/v2.0.1/EzReson_v2.0.1_release.zip
NOTE: For a tutorial in Chinese, a series of articles with examples are available at the WeChat public account, "Quantum Chemistry", following the link:
This WeChat page (in Chinese as well) at "KouShare" provides an overview of the applications of EzReson to diverse chemical problems:
https://mp.weixin.qq.com/s/2JKLrs0UF92Tktst9zWIg
What's new in the latest 2.0.1 version?

Enumeration of all possible Clar structures (with a maximum number of sextets) and all possible Clar resonators (with a variable number of sextets) for pi conjugated compounds

Resonance analysis based on the Clar resonators for pi conjugated compounds

Energy partitioning of the oneelectron energy of each Clar resonator into the contributions from Clar sextets and C=C bonds
How to cite
If you havee used EzReson in your research papers or presentations, it is obligatory to cite the following works:

Yang Wang. A Reliable and Efficient Resonance Theory Based on Analysis of DFT Wave Functions. Phys. Chem. Chem. Phys. 2021, 23, 23312348.

Yang Wang. Superposition of Waves or Densities: Which is the Nature of Chemical Resonance? J. Comput. Chem. 2021, 42, 412417.

Yang Wang. Response to comment on “Superposition of waves or densities: Which is the nature of chemical resonance?” J. Comput. Chem. 2021, 42, 13411343.
If you have used EzReson to generate or analyze Kekulé structures of pi conjugated compounds, it is required to cite the following paper:
 Yang Wang. Extension and Quantification of Fries Rule and Its Connection to Aromaticity: LargeScale Validation by WaveFunctionBased Resonance Analysis. J. Chem. Inf. Model. 2021, in press, (DOI: 10.1021/acs.jcim.1c00735)
If you have used EzReson to generate or analyze Clar structures/resonators of pi conjugated compounds, it is required to cite the following paper:
 Yang Wang. Quantitative Resonance Theory Based on the Clar Sextet Model. J. Phys. Chem. A 2022, 126, 1, 164–176.
Copyright and license
The Author of the EzReson software is Yang Wang ([email protected]; orcid.org/0000000325402199). The EzReson program is released under GNU General Public License v3 (GPLv3).
Disclaimer
The EzReson software is provided as it is, with no warranties. The Author shall not be liable for any use derived from it. Feedbacks and bug reports are always welcome ([email protected]). However, it is kindly reminded that the Author does not take on the responsibility of providing technical support.
How to install
Requirements
 python >= 3.6
 numpy >= 1.18.0
 scipy >= 1.5.1
Installation
 Place the folder of the EzReson package to any location as you like, which is referred to as the source directory hereafter.
 Open with a text editor the script file "ezreson" in the source directory
and set the
EZRESON_DIR
variable as the path of the source directory.  Add the source directory to the global environment variable
PATH
in e.g., ".bash_profile" or ".bashrc" under your HOME directory.
Then, you are ready to go.
How to use
Gaussian calculations
 In the Gaussian input file (e.g., abc.gjf), add in the route section the
keywords
fchk=All Pop=NBO6Read
, and at the end of the file add$NBO NOBOND AONAO=W $END
. In this way, the checkpoint file "Test.FChk" and the NBO matrix file "abc.33" will be generated. Then, rename "Test.FChk" to "abc.fchk". The Gaussian output file should have the extension name of ".out" (If necessary, "abc.log" ought to be renamed as "abc.out").
NOTES:
 It is strongly recommended to use the NBO program later than version 5.0. The free version of NBO 3.1 implemented in Gaussian package would be problematic and give unreliable results.
 Do not use
fchk=All
to generate the checkpoint file if there are two such jobs running at the same working directory at the same time. Otherwise, the two jobs will write the same "Test.FChk" file. Instead, add the%chk=abc.chk
line to obtain the checkpoint file "abc.chk". Then, use Gaussian's formchk utility to convert "abc.chk" to "abc.fchk".
 Make sure that the following four files, as the inputs for EzReson, are in the same working directory:
 abc.gjf
 abc.out
 abc.fchk
 abc.33

Change to the working directory and prepare the input file for EzReson (vide infra), say, "abc_wfrt.in"

Simply execute the following command:
ezreson abc_wfrt.in > abc_wfrt.out
As you see, you will find the result of resonance analysis by EzReson in file "abc_wfrt.out".
NOTE: If the analysis is conducted within the framework of the simple Hückel molecular orbital (HMO) theory, only an xyz or gjf file is required to feed EzReson.
A typical wavefunctionbased resonance theory (WFRT) analysis
Let us take benzene as an example. Suppose that after the DFT calculations of benzene you have obtained the four output files, Ph.gjf, Ph.out, Ph.fchk and Ph.33.
1. Perform an LMO calculation to obtain PipekMezey localized MOs
Prepare an input file, named "benzene_lmo.in" as:
File = Ph
Job = LMO
NOTE: The letters in EzReon's input file are caseinsensitive.
Then, use the following command to run the LMO job:
ezreson benzene_lmo.in > benzene_lmo.out
After finishing the LMO calculation correctly, the following output files will be generated:
 Ph_CNAOLMO.dat
 Ph_ELMO.dat
 Ph_LMO.fchk
2. Identify the LMOs associated with the resonance subsystem
In this particular case of benezene, we are to find the occupied LMOs
corresponding to the pi conjugated system.
To this end, open file "Ph_LMO.fchk" with visualization software like JMol or
Gabedit. For JMol, after opening Ph_LMO.fchk, type in the script console:
isosurface mo 21
and the LMO#21 (which is the HOMO) will be displayed. You
will see that this LMO belongs to the piresonance system and is thus selected.
Then, keep on inspecting lower LMOs, #20, #19, ..., until all LMOs belonging to
the resonance subsystem have been chosen. Since in this particular case there
are 6 electrons in the resonance subsystem, you only need to identify 3 LMOs.
As for benzene, the resonating LMOs are identified as 19, 20 and 21.
NOTE: Other visualization softwares may not support reading *.fchk for visualization of orbitals. But you can always convert *.fchk file to *.cube file by the cubegen utility in the Gaussian suite. Then, a great variety of visualization tools are able to read the *.cube file for visualization of molecular orbitals.
3. Perform the WRFT analysis
We have determined that the LMOs for the resonance subsystem of benzene are #19, 20 and 21. We also see that the involved atoms are the carbon atoms, whose indices are 1, 2, 3, 4, 5 and 6 as indicated in "Ph.gjf".
Accordingly, we prepare an input file, named "benzene_wfrt.in" as:
File = Ph
Job = WFRT
LMOs = 19 20 21
Atoms = 1 2 3 4 5 6
Then, use the following command to run the WFRT job:
ezreson benzene_wfrt.in > benzene_wfrt.out
NOTE: i) For the indices of atoms, the order matters in order to apply Rumer's rule for determination of linearly independent set of Lewis structures. For monocyclic systems, the ordered atoms should form a circle. For other systems, the choice is somewhat arbitrary, but it is recommended that the atoms be disposed to form a circle as much as possible.
ii) For versions 1.1.4+, the LMOs
and Atoms
can be defined in a more compact
way by using the MATLAB colon syntax:
LMOs = 19:21
Atoms = 1:6
iii) In this case of benzene (and pi conjugated systems in general), the pi LMOs
can be automatically chosen by EzReson with the PI
keyword (caseinsensitive):
LMOs = PI
Note that this is only valid for allcarbon conjugated systems, and not for heteroatomic systems yet.
Fragment of the output:
Performing Wave Function based Resonance Theory (WFRT) calculations ...
Reading LEWIS_6c_6e.dat for Lewis structures ...
215 Lewis structures in total to be considered
... ...
... ...
Solving linear equation for WFRT coefficients ...
40 Lewis structures violating Rumer's rule excluded
0 Lewis structures of small projections ( < 0.00E+00 ) excluded
Using totally 175 Lewis structures to expand the actual wave function
Normalizing the wave function by a factor of 0.9951914286 ...

No. Projection Coefficient RE Mulli. Bickel. RosSc. Lowdin PWSO Lewis structure

1 0.559795 0.3333335 175.96 18.75% 22.85% 9.60% 14.59% 19.80% 12 34 56
2 0.559795 0.3333335 175.96 18.75% 22.85% 9.60% 14.59% 19.80% 23 45 16
15 0.221154 0.2222227 347.79 4.94% 2.70% 4.27% 2.77% 2.44% 4: 23 56
16 0.221154 0.2222227 347.79 4.94% 2.70% 4.27% 2.77% 2.44% 1: 23 56
17 0.221154 0.2222219 347.80 4.94% 2.70% 4.27% 2.77% 2.44% 6: 12 45
18 0.221154 0.2222219 347.80 4.94% 2.70% 4.27% 2.77% 2.44% 3: 12 45
19 0.221154 0.2222219 347.80 4.94% 2.70% 4.27% 2.77% 2.44% 2: 34 16
20 0.221154 0.2222219 347.80 4.94% 2.70% 4.27% 2.77% 2.44% 5: 34 16
... ...
... ...

Reproducibility = 99.519143%
WFRT analysis using Lewis structures with maximum number of lone pairs
Use the option MAXNLP = <num>
in the input file to set the maximum number of
lone pairs allowed in the Lewis structures
NOTE: <num>
can be 0, 1, 2, ..., or inf
WFRT analysis using projection cutoff to truncate the set of Lewis structures
For instance, if you want to employ only the Lewis structures whose projection to the actual wave function is less than 0.1, just add the following line in the input file:
ProjCut = 0.1
WFRT analysis using specified Lewis structures
You can perform the WFRT analysis using only a userdefined Lewis structures,
which can be specified in the input file by using the Lewis
keyword. There are
two ways of specifying Lewis structures.

Specification with indices of Lewis structures:
Here is an example:
Lewis = 1 2 9 10
It means that only the 1st, 2nd, 9th and 10th Lewis structures are to be used to expand the wave function. Note that the indices follow the descending order of the projection of the Lewis structures onto the actual wave function.

Specification by writing explicitly the Lewis structures
An example:
Lewis = 12/34/56 1:/23/5: 2:/4:/6:
In this example, we have defined three Lewis structures. The first one is a Kekulé structure, where

denotes the bond connecting two atoms. In the second Lewis structure, each of atoms 1 and 5 has a lone pair and atoms 2 and 3 forms a bond. The last Lewis structure contains three lone pairs located, respectively, at atoms 2, 4 and 6.
WFRT analysis using all Kekulé structures
To perform a WFRT analysis using only all possible Kekulé structures, set the
Kekule
parameter to be TRUE:
Kekule = TRUE
NOTE:
i) The two options Kekule
and Lewis
are repulsive to each other,
i.e., they cannot be present at the same time in the input file.
ii) Before performing the WFRT analysis, Kekulé structures are first enumerated
and stored in a file with *.kek extension. You can also generate the Kekulé
structures without performing WFRT or PROJ job, by invoking the ENUM
job type:
Job = ENUM
WFRT analysis using all Clar resonators
To perform a WFRT analysis using all possible Clar resonators, set the Clar
parameter to be TRUE:
Clar = TRUE
If you want to perform a WFRT analysis using only the Clar structures (i.e.,
with the maximum possible number of Clar sextets), set the ClarMax
parameter to be TRUE:
ClarMax = TRUE
NOTE:
i) The options Clar
, ClarMax
and Kekule
are mutually repulsive to each other, i.e., they cannot be present at the same time in the input file.
ii) The Job types, PROJ
and ENUM
, and the Huckel
option are also valid for Clar and ClarMax analysis. But the DMRT
analysis is not supported for Clartype resonance structures.
WFRT analysis in the framework of simple Hueckel molecular orbital (HMO) theory
The WFRT analysis can be performed in the simple HMO framework. This is done by
turning on the Huckel
option in the input file, and there is no need to set
the LMOs.
The simplest example:
File = Ph
Job = WFRT
Atoms = 1 2 3 4 5 6
Huckel = TRUE
Furthermore, the use of Huckel
can also be accompanied with other keywords,
including MaxNLP
, ProjCut
, Kekule
and Lewis
, for the respective
purposes.
A typical densitymatrixbased resonance theory (DMRT) analysis
A typical example for a DMRT analysis is as follows:
file = Ph
job = DMRT
LMOs = 19 20 21
atoms = 1 2 3 4 5 6
The keywords such as Lewis
and Kekule
are also valid for DMRT analysis.
NOTES:

We strongly discourage the application of the DMRT, for this theory has been proved to have fundamental and intrinsic inadequacies.

One may still use the projections of the density matrices associated with the Lewis strutures onto the actual density matrix of the resonance subsystem, as they provide a way to measure how a given Lewis struture ressembles the actual electronic structure. For this purpose, a
PROJ
job should be carried out (see below).
DMRT analysis in the HMO framework
Just switch on the Huckel
option to perform a DMRT analysis in the HMO
framework.
Projection calculations
The calculation of projections is much cheaper than a full expansion calculation, as the former scales as O(N) while the latter scales from O(N^3) to O(N^2). Projecitons can also provide useful information. The wave function or density matrix projection of a Lewis structure may be regarded as the resemblence of that Lewis structure to the actual molecular structure.
To launch a projection job, just indicate it in the Job
parameter:
File = Ph
Job = PROJ
LMOs = 19 20 21
Atoms = 1 2 3 4 5 6
By running this job, wave function projections, density matrix projections and the energies of Lewis structures will be obtained in the same output.
Other options, like Lewis
, Huckel
can be used for a PROJ
job.
Cautions
It is recommended to check the phase matching between the reduced atomic orbitals (RAOs) if the analysis is not on the basis of HMO theory and the resonance system is not a planar pi conjugated system.
Usually, the phases of the RAOs are automatically determined by EzReson, which works quite successfully for most of the cases. However, somehow it might fail occasionally. A known example is given by the 3c4e pi conjugated system in the ozone molecule. In the HOMO of O3, the two terminal O atoms have an offphase combination, which is rather exceptional.
To verify that all RAO phases automatically determined by EzReson are correctly
matched, one should active the writeRAOs
option in the input file to output
the RAOs:
file = Ph
job = PROJ
LMOs = 19 20 21
atoms = 1 2 3 4 5 6
Lewis = 1
writeRAOs = TRUE
This is a simple projection job that only calculate the projection of one Lewis
structure (so it is very fast). By running this job, a file named "Ph_RAO.fchk"
is generated. By opening this fchk with softwares like JMol, one can visulized
the RAOs atom by atom to check if all RAOs are in phase.
What to do if all RAOs are not in phase?

Set the
flipRAOs
to zero and rerun the job:writeRAOs = TRUE flipRAOs = 0
and obtain a new fchk file of RAOs, "Ph_RAO.fchk".

Visualize the new set of RAOs and pick out all RAOs with opposite phases (just record the corresponding indices of atoms, say, 2, 3, 6)

Set the
flipRAOs
to the indices of atoms whose RAO needs to be flipped:writeRAOs = TRUE flipRAOs = 2, 3, 6

Rerun the job and visualize "Ph_RAO.fchk" again to finally verify the phase matching of the RAOs.
Limitations
The current WFRT scheme have some limitations in its applicability. In the near future, we will be endeavoring to extend and generalize the approach for a wider range of applications.
The following aspects are to be considered in our future work:

Openshell systems

Going beyond the effective minimal atomic basis approximation, which will improve the reproducibilty of wave function

Using the varying RAOs that depend on each specific Lewis structures and allowing to use multiple RAOs for one atom, which will make the method able to deal with systems like SO_4^{2}, PF5, B2H6, etc.

Going beyond DFT/HF wave functions to enable the expansion of multireference wave functions

Excited states