Getting started¶

In this section we will go through the process of performing all the available calculations using ESIpy’s utilities based on the provided examples. There are two main ways to run ESIpy: by creating a PySCF input, running the single-point calculation and running ESIpy afterwards, or by creating an input from a FCHK file. The second option may not contain all the features, they will be added in future releases.

ESIpy from PySCF¶

ESIpy works on the object ESI, which will contain all the information required for the calculation. It is recommended to initialize the object with all the data, rather than adding it once the initialization process is finished. To run the code from the terminal, generate the Python script or adapt those of this webpage and run it as python3 code.py. To save the output as a file, use python3 code.py > code.esi. Even though the .esi extension is not mandatory, we recommend using it. The test files have been generated from benzene and naphthalene at singlet (and triplet) level of theory, and are just set to showcase the capabilities of the code.

The simplest form of input following a usual PySCF calculation

from pyscf import gto, dft
import esipy

mol = gto.Mole()
mol.atom = '''
6        0.000000000      0.000000000      1.393096000
6        0.000000000      1.206457000      0.696548000
6        0.000000000      1.206457000     -0.696548000
6        0.000000000      0.000000000     -1.393096000
6        0.000000000     -1.206457000     -0.696548000
6        0.000000000     -1.206457000      0.696548000
1        0.000000000      0.000000000      2.483127000
1        0.000000000      2.150450000      1.241569000
1        0.000000000      2.150450000     -1.241569000
1        0.000000000      0.000000000     -2.483127000
1        0.000000000     -2.150450000     -1.241569000
1        0.000000000     -2.150450000      1.241569000
'''
mol.basis = 'sto-3g'
mol.spin = 0
mol.charge = 0
mol.symmetry = True
mol.verbose = 9
mol.build()

mf = dft.KS(mol)
mf.kernel()

ring = [1, 2, 3, 4, 5, 6]
arom = esipy.ESI(mol=mol, mf=mf, rings=ring, partition="nao")
arom.print()

By providing the mol and mf objects, ESIpy generates the AOMs in the desired partition and computes the indices following the ring connectivity in the list rings. The program works similarly through unrestricted wavefunctions, the output of which provides the indices split into orbital contributions.

Note

In the following, we will only consider the ESIpy part of the code.

Dealing with AOMs¶

In order to avoid the single-point calculation, the save attribute will save the AOMs and a dictionary containing information about the molecule and calculation into a binary file in disk. There will be a directory created with the name of the calculation and the _esipy/ extension, and all the .aoms and .molinfo files will be saved there. It should contain only the title of the calculation, and ESIpy will add the used partition scheme and the extension (.aoms and .molinfo). Hereafter, these will be accessible at any time. It is also recommended to use a for-loop scheme for all the partitions, as the computational time to generate the matrices is minimal and independent to the chosen scheme.

ring = [1, 2, 3, 4, 5, 6]
name = "benzene"
for part in ["mulliken", "lowdin", "meta-lowdin", "nao", "iao"]:
    arom = esipy.ESI(mol=mol, mf=mf, rings=ring, partition=part, save=name)
    arom.print()

Additionally, one can generate a directory containing the AOMs in AIMAll format. These files are readable from ESIpy, but also from Eduard Matito’s ESI-3D code. These are written through the method writeaoms()

arom = esipy.ESI(mol=mol, mf=mf, rings=[1,2,3,4,5,6], partition="nao")
arom.writeaoms("benzene_nao")

and read through the method readaoms() by previously specifying the read=True attribute and using another generated “molinfo” dictionary.

arom = esipy.ESI(rings=[1,2,3,4,5,6], partition="nao", read=True, molinfo="benzene_nao.molinfo", aoms="benzene_nao.aoms")
arom.readaoms()
arom.print()

Correlated wavefunctions¶

For natural orbitals wavefunctions, an additional diagonalization of the first-order reduced density matrix (1-RDM) is carried out, the computational time of which is also very low. The single-determinant (RHF or UHF) object has to be provided through the myhf attribute. Population analyses use both Fulton’s approach and the 2-RDM approximation in terms of natural occupations, but only Fulton’s approach is used for the aromaticity calculations.

from pyscf import gto, scf, ci, cc, mp, mcscf
import esipy

mol = gto.Mole()
mol.atom = '''
6        0.000000000      0.000000000      1.393096000
6        0.000000000      1.206457000      0.696548000
6        0.000000000      1.206457000     -0.696548000
6        0.000000000      0.000000000     -1.393096000
6        0.000000000     -1.206457000     -0.696548000
6        0.000000000     -1.206457000      0.696548000
1        0.000000000      0.000000000      2.483127000
1        0.000000000      2.150450000      1.241569000
1        0.000000000      2.150450000     -1.241569000
1        0.000000000      0.000000000     -2.483127000
1        0.000000000     -2.150450000     -1.241569000
1        0.000000000     -2.150450000      1.241569000
'''
mol.basis = 'sto-3g'
mol.spin = 0
mol.charge = 0
mol.symmetry = True
mol.verbose = 0
mol.max_memory = 4000
mol.build()

mf = scf.RHF(mol).run()

print("Running CCSD calculation...")
mf1 = cc.CCSD(mf).run()
print("Running CISD calculation...")
mf2 = ci.CISD(mf).run()
print("Running CASSCF calculation...")
mf3 = mcscf.CASSCF(mf, 6, 6).run()
print("Running MP2 calculation...")
mf4 = mp.MP2(mf).run()
ring = [1, 2, 3, 4, 5, 6]

for part in ["mulliken", "lowdin", "meta-lowdin", "nao", "iao"]:
    for method in [mf1, mf2, mf3, mf4]:
        arom = esipy.ESI(mol=mol, mf=method, myhf=mf, rings=ring, partition=part)
        arom.print()

Note

The IAOs expand the occupied orbitals in the same rank as the minimal basis, but the role of valence orbitals is important for the calculation. Therefore, the transformation matrix is computed through the RHF object, thus making the myhf attribute needed for these calculations. However, it is recommended to use other robust schemes for multi-determinant wave functions.

ESIpy from a FCHK file¶

ESIpy can also read wavefunctions from a FCHK file generated by Gaussian. This is done by generating a simple input. As an example, the simplest form of input would be:

$READFCHK
benzene.fchk
$RING
1 2 3 4 5 6
$PARTITION
META_LOWDIN
NAO
IAO

All the keywords available are explained in the ESIpy FCHK input reference section.