Tutorial 1: Introduction to Antimatter Quantum Chemistry

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
from antinature.core.basis import MixedMatterBasis, GaussianBasisFunction
from antinature.core.hamiltonian import AntinatureHamiltonian
from antinature.core.integral_engine import AntinatureIntegralEngine
from antinature.specialized import PositroniumSCF
from antinature.core.molecular_data import MolecularData

# Create a minimal basis set manually
basis = MixedMatterBasis()
center = np.array([0.0, 0.0, 0.0])  # Origin

# Create s-type electron basis functions
e_basis = []
e_basis.append(GaussianBasisFunction(center, 1.0, (0, 0, 0)))
e_basis.append(GaussianBasisFunction(center, 0.5, (0, 0, 0)))

# Create s-type positron basis functions (slightly more diffuse than electrons)
p_basis = []
p_basis.append(GaussianBasisFunction(center, 0.8, (0, 0, 0)))
p_basis.append(GaussianBasisFunction(center, 0.4, (0, 0, 0)))

# Set up the basis
basis.electron_basis.basis_functions = e_basis
basis.positron_basis.basis_functions = p_basis
basis.n_electron_basis = len(e_basis)
basis.n_positron_basis = len(p_basis)
basis.n_total_basis = basis.n_electron_basis + basis.n_positron_basis

print(f"Created a basis set with {basis.n_electron_basis} electron and {basis.n_positron_basis} positron functions")
print("Electron exponents: 1.0, 0.5")
print("Positron exponents: 0.8, 0.4")
print("(Positron basis functions are more diffuse to represent the extended nature of positronium)")

# Define the positronium system
ps_data = MolecularData.positronium()
print("Positronium system created:")
print(f"- Number of electrons: {ps_data.n_electrons}")
print(f"- Number of positrons: {ps_data.n_positrons}")
print(f"- System charge: {ps_data.charge}")
print(f"- Description: {ps_data.description}")

# Set up integral engine for the basis
integral_engine = AntinatureIntegralEngine(use_analytical=True)
basis.set_integral_engine(integral_engine)

# Create a simplified test Hamiltonian with predefined values
n_e_basis = basis.n_electron_basis
n_p_basis = basis.n_positron_basis
n_total = n_e_basis + n_p_basis

# Create overlap matrix (S)
S = np.zeros((n_total, n_total))
S[:n_e_basis, :n_e_basis] = np.array([[1.0, 0.8], [0.8, 1.0]])  # Electron block
S[n_e_basis:, n_e_basis:] = np.array([[1.0, 0.7], [0.7, 1.0]])  # Positron block

# Create core Hamiltonian matrices (H_core)
H_core_e = np.array([[-1.0, -0.5], [-0.5, -0.8]])  # For electrons
H_core_p = np.array([[-1.0, -0.4], [-0.4, -0.7]])  # For positrons

# Create electron-positron attraction tensor (ERI_ep)
ERI_ep = np.zeros((n_e_basis, n_e_basis, n_p_basis, n_p_basis))
for i in range(n_e_basis):
    for j in range(n_e_basis):
        for k in range(n_p_basis):
            for l in range(n_p_basis):
                # Simple exponential decay model for attraction
                ERI_ep[i, j, k, l] = -0.5 * np.exp(-((i - k) ** 2 + (j - l) ** 2))

# Create Hamiltonian dictionary
hamiltonian_matrices = {
    'overlap': S,
    'H_core_electron': H_core_e,
    'H_core_positron': H_core_p,
    'electron_positron_attraction': ERI_ep,
}

print("Hamiltonian components created:")
print("- Overlap matrix (S): 4×4 matrix")
print("- Electron core Hamiltonian: 2×2 matrix")
print("- Positron core Hamiltonian: 2×2 matrix")
print("- Electron-positron attraction: 2×2×2×2 tensor")

# Create the SCF solver for para-positronium
ps_scf = PositroniumSCF(
    hamiltonian=hamiltonian_matrices,
    basis_set=basis,
    molecular_data=ps_data,
    positronium_state='para',  # para-positronium (singlet state)
    include_qed_corrections=True,  # Include quantum electrodynamics corrections
    convergence_threshold=1e-5,
    damping_factor=0.5,
    max_iterations=100
)

# Solve the SCF equations
results = ps_scf.solve_scf()

# Print the results
print(f"Energy: {results['energy']:.6f} Hartree")
print(f"Theoretical energy: -0.25 Hartree")
print(f"Iterations: {results.get('iterations', 'N/A')}")
print(f"Converged: {results.get('converged', 'N/A')}")

if 'exact_solution' in results and results['exact_solution']:
    print("Used exact analytical solution for positronium")

# Calculate annihilation properties
try:
    annihilation_data = ps_scf.calculate_annihilation_rate()
    print(f"Annihilation rate: {annihilation_data['rate']:.4e} s^-1")
    print(f"Lifetime: {annihilation_data['lifetime']:.4e} s")
    print("(Para-positronium has a lifetime of ~125 picoseconds in vacuum)")
except Exception as e:
    print(f"Could not calculate annihilation properties: {str(e)}")
    print("This is expected with our simplified Hamiltonian")
    print("Theoretical lifetime of para-positronium: ~125 picoseconds")

# Create ortho-positronium SCF solver
ortho_ps = PositroniumSCF(
    hamiltonian=hamiltonian_matrices,
    basis_set=basis,
    molecular_data=ps_data,
    positronium_state='ortho',  # ortho-positronium (triplet state)
    include_qed_corrections=True,
    convergence_threshold=1e-5,
    damping_factor=0.5,
    max_iterations=100
)

# Solve SCF for ortho-positronium
ortho_results = ortho_ps.solve_scf()

# Print comparison
print("\n== Comparison of Para vs. Ortho Positronium ==")
print(f"Para-positronium energy: {results['energy']:.6f} Hartree")
print(f"Ortho-positronium energy: {ortho_results['energy']:.6f} Hartree")

# Calculate annihilation rate for ortho-positronium
try:
    ortho_annihilation = ortho_ps.calculate_annihilation_rate()
    print("\nAnnihilation properties comparison:")
    print(f"Para-positronium lifetime: {annihilation_data['lifetime']:.4e} s")
    print(f"Ortho-positronium lifetime: {ortho_annihilation['lifetime']:.4e} s")
    print(f"Lifetime ratio (ortho/para): {ortho_annihilation['lifetime']/annihilation_data['lifetime']:.2f}")
except Exception as e:
    print("\nCould not calculate ortho-positronium annihilation properties")
    print("Theoretical lifetimes:")
    print("- Para-positronium: ~125 picoseconds")
    print("- Ortho-positronium: ~142 nanoseconds")
    print("- Ratio: ~1100")

try:
    ps_scf.visualize_orbitals(
        grid_dims=(20, 20, 20),
        limits=(-5.0, 5.0),
        save_path="positronium_orbital.png"
    )
    print("Visualization saved to 'positronium_orbital.png'")
except Exception as e:
    print(f"Visualization error: {str(e)}")
    print("Positronium wavefunctions typically have spherical symmetry around the center of mass")