Quick Start

This notebook provides a quick introduction to using MolGrid for molecular grid generation and integration. We’ll walk through the basic workflow of creating molecules, generating grids, and calculating electron density with PySCF.

Creating a Molecule

First, let’s create a simple water molecule using the Molecule class. The water molecule has an oxygen atom at the origin and two hydrogen atoms at specific positions.

from molgrid import Atom, Molecule

# Create atoms with their coordinates (in angstroms)
oxygen = Atom('O', [0.0, 0.0, 0.0])
hydrogen1 = Atom('H', [0.0, -0.757, 0.587])
hydrogen2 = Atom('H', [0.0, 0.757, 0.587])

# Create molecule from the atoms
water = Molecule([oxygen, hydrogen1, hydrogen2])
print(f"Created molecule: {water}")
print(f"Number of atoms: {len(water)}")

# Let's also print the atom information
print("\nAtom details:")
for i, atom in enumerate(water):
    print(f"  Atom {i}: {atom.symbol} at {atom.coordinate}")

Created molecule: <molgrid.molecule.Molecule object at 0x11a84f390>
Number of atoms: 3

Atom details:
  Atom 0: O at [0. 0. 0.]
  Atom 1: H at [ 0.    -0.757  0.587]
  Atom 2: H at [0.    0.757 0.587]

Creating an Atomic Grid

Now let’s create an atomic grid for the oxygen atom. Atomic grids are constructed using radial shells and angular points. The radial shells determine how many concentric spheres we have, and the angular points determine how many points on each sphere.

from molgrid import AtomicGrid

# Create atomic grid for oxygen
# nshells=10: 10 radial shells
# nangpts=110: 110 angular points per shell (using Lebedev quadrature)
oxygen_grid = AtomicGrid(oxygen, nshells=75, nangpts=302)

print(f"Oxygen atomic grid:")
print(f"  Number of points: {len(oxygen_grid)}")
print(f"  Coordinates shape: {oxygen_grid.coords.shape}")
print(f"  Weights shape: {oxygen_grid.weights.shape}")
print(f"  Total weight: {oxygen_grid.weights.sum():.6f}")

# Let's look at the first few points
print("\nFirst 5 grid points and weights:")
for i in range(5):
    print(f"  Point {i}: {oxygen_grid.coords[i]}, weight: {oxygen_grid.weights[i]:.6f}")

Oxygen atomic grid:
  Number of points: 22650
  Coordinates shape: (22650, 3)
  Weights shape: (22650,)
  Total weight: 59199645676.309372

First 5 grid points and weights:
  Point 0: [1326.73399783    0.            0.        ], weight: 50173399.244666
  Point 1: [-1326.73399783     0.             0.        ], weight: 50173399.244666
  Point 2: [   0.         1326.73399783    0.        ], weight: 50173399.244666
  Point 3: [    0.         -1326.73399783     0.        ], weight: 50173399.244666
  Point 4: [   0.            0.         1326.73399783], weight: 50173399.244666

Creating a Molecular Grid

A molecular grid combines atomic grids from all atoms in the molecule. We use Becke partitioning to assign weights to each grid point based on its proximity to different atoms. This prevents double-counting of electron density.

from molgrid import MolecularGrid

# Create molecular grid with Becke partition
# partition_method='becke': Use Becke partitioning for grid partition
mol_grid = MolecularGrid(water, nshells=75, nangpts=302, partition_method='becke')

print(f"Molecular grid:")
print(f"  Number of points: {len(mol_grid)}")
print(f"  Coordinates shape: {mol_grid.coords.shape}")
print(f"  Weights shape: {mol_grid.weights.shape}")
print(f"  Total weight: {mol_grid.weights.sum():.6f}")

# Compare with unpartitioned grid
mol_grid_no_prune = MolecularGrid(water, nshells=75, nangpts=302, partition_method=None)
print(f"\nWithout partitioning: {len(mol_grid_no_prune)} points")
print(f"  Total weight: {mol_grid_no_prune.weights.sum():.6f}")

Molecular grid:
  Number of points: 67950
  Coordinates shape: (67950, 3)
  Weights shape: (67950,)
  Total weight: 28788884690.918945

Without partitioning: 67950 points
  Total weight: 82701356864.936829
  Total weight: 82701356864.936829

Iterating Through Atomic Grids

You can iterate through the atomic grids that make up the molecular grid. This is useful for analyzing contributions from individual atoms.

print("Atomic grids in pruned molecular grid:")
total_points = 0
for i, atomic_grid in enumerate(mol_grid):
    points = len(atomic_grid)
    total_points += points
    print(f"  Atom {i} ({atomic_grid.atom.symbol}): {points} points, weight sum: {atomic_grid.weights.sum():.6f}")

print(f"\nTotal points from atomic grids: {total_points}")
print(f"Matches molecular grid total: {total_points == len(mol_grid)}")

Atomic grids in pruned molecular grid:
  Atom 0 (O): 22650 points, weight sum: 21244363291.203724
  Atom 1 (H): 22650 points, weight sum: 3772260699.857614
  Atom 2 (H): 22650 points, weight sum: 3772260699.857613

Total points from atomic grids: 67950
Matches molecular grid total: True

Calculating Electron Density

Now let’s calculate the electron density using PySCF. We’ll use the molecular grid we created to integrate the electron density and get the total number of electrons in the molecule.

import numpy as np
from pyscf import gto, dft
from pyscf.dft import numint

# Create PySCF molecule (same geometry as our MolGrid molecule)
mol = gto.M(atom='O 0 0 0; H 0 -0.757 0.587; H 0 0.757 0.587', basis='6-31g', verbose=0)

# Run DFT calculation to get the electron density
mf = dft.RKS(mol)
mf.kernel()
rdm1 = mf.make_rdm1()  # Get the 1-electron reduced density matrix

# Evaluate atomic orbitals on our grid
ao = numint.eval_ao(mol, mol_grid.coords, deriv=0)

# Calculate electron density at each grid point
rho = numint.eval_rho(mol, ao, rdm1, xctype='LDA')

# Integrate electron density over the entire grid
number_of_electrons = np.sum(rho * mol_grid.weights)

print(f"Electron density calculation results:")
print(f"  Calculated number of electrons: {number_of_electrons:.6f}")
print(f"  Expected number of electrons: 10.000000")
print(f"  Error: {abs(number_of_electrons - 10.0):.6f}")

# Let's also calculate electron density for each atomic grid
print("\nElectron distribution by atom:")
total_electrons = 0
for i, atomic_grid in enumerate(mol_grid):
    # Evaluate AO on atomic grid
    ao_atom = numint.eval_ao(mol, atomic_grid.coords, deriv=0)
    # Calculate density on atomic grid
    rho_atom = numint.eval_rho(mol, ao_atom, rdm1, xctype='LDA')
    # Integrate electron density on atomic grid
    atom_electrons = np.sum(rho_atom * atomic_grid.weights)
    total_electrons += atom_electrons
    print(f"  Atom {i} ({atomic_grid.atom.symbol}): {atom_electrons:.6f} electrons")

print(f"\nTotal electrons from atomic grids: {total_electrons:.6f}")
print(f"Matches molecular total: {abs(total_electrons - number_of_electrons) < 1e-10}")

Electron density calculation results:
  Calculated number of electrons: 10.000001
  Expected number of electrons: 10.000000
  Error: 0.000001

Electron distribution by atom:
  Atom 0 (O): 6.820235 electrons
  Atom 1 (H): 1.589883 electrons
  Atom 2 (H): 1.589883 electrons

Total electrons from atomic grids: 10.000001
Matches molecular total: True