Calculating the Mulliken population

Calculating the Mulliken population#

In Cedillo2012, we find a Mulliken charge of 0.20 on the carbon atom in CO. In this notebook, we will reproduce this result.

# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')

import gqcpy
import numpy as np

Let’s start by creating the molecule, and quantizing the Hamiltonian in the associated restricted spin-orbital basis of the AOs.

Restricted Mulliken populations#

molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/CO_mulliken.xyz" , 0)  # '0': Create a neutral molecule.
N = molecule.numberOfElectrons()

spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
S = spin_orbital_basis.quantize(gqcpy.OverlapOperator())

fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
hamiltonian = spin_orbital_basis.quantize(fq_hamiltonian)

Afterwards, we should solve the RHF SCF equations.

environment = gqcpy.RHFSCFEnvironment_d.WithCoreGuess(N, hamiltonian, S)
solver = gqcpy.RHFSCFSolver_d.DIIS()
objective = gqcpy.DiagonalRHFFockMatrixObjective_d(hamiltonian)  # This objective makes sure that the RHF parameters represent the canonical RHF MOs.

rhf_parameters = gqcpy.RHF_d.optimize(objective, solver, environment).groundStateParameters()
D = rhf_parameters.calculateScalarBasis1DM()  # The AO density matrix.

The Mulliken population is effectively the expectation value of the Mulliken operator and a density matrix, so we’ll have to set up the Mulliken operator.

Since we’re interested in determining the Mulliken charge on the carbon atom, we’ll prepare a Mulliken partitioning of the spin-orbital basis with basis functions centered on the carbon atom.

mulliken_domain = spin_orbital_basis.mullikenDomain(lambda shell: shell.nucleus().element() == "C")

Let’s also confirm that the supplied selector lambda function has selected the correct indices for the basis functions centered on the carbon atom.

print(mulliken_domain.domainIndices())

[0, 1, 2, 3, 4]

Theoretically, the Mulliken operator is the Mulliken-partitioned number/overlap operator. Since GQCP supports partitioning of any one-electron operator, we can use the .partitioned() API.

P = mulliken_domain.projectionMatrix(spin_orbital_basis.expansion())
mulliken_operator = S.partitioned(P)

What remains is to calculate the expectation value of the Mulliken-partitioned number operator, which yields the Mulliken population. The Mulliken charge can then be calculated straightforwardly.

C_population = mulliken_operator.calculateExpectationValue(D)[0]
print("Mulliken population: ", C_population)

Mulliken population:  5.801066370434745

C_charge = 6 - C_population
print("Mulliken charge: ", C_charge)

Mulliken charge:  0.1989336295652553

Unrestricted Mulliken populations#

GQCP also offers the ability to calculate alpha- and beta- Mulliken populations. Since the previous example would be trivial, let’s use CO+ instead of CO.

molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/CO_mulliken.xyz" , +1)  # '+1': Create CO+.
N_alpha = molecule.numberOfElectronPairs() + 1
N_beta = molecule.numberOfElectronPairs()

u_spin_orbital_basis = gqcpy.USpinOrbitalBasis_d(molecule, "STO-3G")
S_unrestricted = u_spin_orbital_basis.quantize(gqcpy.OverlapOperator())

fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
hamiltonian = u_spin_orbital_basis.quantize(fq_hamiltonian)  # We perform UHF calculations with the Hamiltonian in the basis of the restricted atomic spin-orbitals. 

We solve the UHF SCF equations.

environment = gqcpy.UHFSCFEnvironment_d.WithCoreGuess(N_alpha, N_beta, hamiltonian, S_unrestricted)
solver = gqcpy.UHFSCFSolver_d.DIIS(threshold=1.0e-6, maximum_number_of_iterations=500)

uhf_parameters = gqcpy.UHF_d.optimize(solver, environment).groundStateParameters()
D_unrestricted = uhf_parameters.calculateScalarBasis1DM()  # The AO density matrix.

We now analogously proceed by setting up the Mulliken partitioning scheme, and finally partition the number/overlap operator according to the obtained partitioning scheme.

u_mulliken_domain = u_spin_orbital_basis.mullikenDomain(lambda shell: shell.nucleus().element() == "C")
u_P = u_mulliken_domain.projectionMatrix(u_spin_orbital_basis.expansion())
u_mulliken_operator = S_unrestricted.partitioned(u_P)

Let’s first check the total Mulliken population and charge.

C_population_total = u_mulliken_operator.calculateExpectationValue(D_unrestricted)[0]
print("Mulliken population: ", C_population_total)
print("Mulliken charge: ", 6 - C_population_total)

Mulliken population:  5.202385763365194
Mulliken charge:  0.7976142366348062

Finally, let’s find out the Mulliken populations for the alpha and beta electrons, separately

C_population_alpha = u_mulliken_operator.alpha.calculateExpectationValue(D_unrestricted.alpha)[0]
C_population_beta = u_mulliken_operator.beta.calculateExpectationValue(D_unrestricted.beta)[0]

print("Mulliken population alpha: ", C_population_alpha)
print("Mulliken population beta: ", C_population_beta)

Mulliken population alpha:  3.5902074357563403
Mulliken population beta:  1.6121783276088537

Generalized Mulliken populations#

GQCP also allows for Mulliken analysis in case of complex orbital bases. The example below applies Mulliken population analysis on a complex GHF calculation.

molecule = gqcpy.Molecule.ReadXYZ("../../gqcp/tests/data/CO_mulliken.xyz" , 0)  # '0': Create a neutral molecule.
N = molecule.numberOfElectrons()

spinor_basis = gqcpy.GSpinorBasis_cd(molecule, "STO-3G")
S_generalized = spinor_basis.quantize(gqcpy.OverlapOperator())

fq_hamiltonian = gqcpy.FQMolecularHamiltonian(molecule)
gsq_hamiltonian = spinor_basis.quantize(fq_hamiltonian) 

We solve the GHF SCF equations.

ghf_environment = gqcpy.GHFSCFEnvironment_cd.WithComplexlyTransformedCoreGuess(N, gsq_hamiltonian, S_generalized)
ghf_solver = gqcpy.GHFSCFSolver_cd.Plain(threshold=1.0e-6, maximum_number_of_iterations=1500)

ghf_parameters = gqcpy.GHF_cd.optimize(ghf_solver, ghf_environment).groundStateParameters()
D_generalized = ghf_parameters.calculateScalarBasis1DM()  # The AO density matrix.