Die Jagd nach dem perfekten Leiter: Ein Blick in die Welt der Hochtemperatur-Supraleiter und warum Quantenchemie hilft

(Ich habe Gemini überzeugt, so eine grobe Abschätzung mit mir zu machen. Gemini hat natürlich immer behauptet, man muss es wie in den Artikeln machen, die dazu veröffentlicht sind)

Stellen Sie sich eine Welt vor, in der elektrische Energie ohne Verluste transportiert wird. Keine Wärmeentwicklung, keine Energieverschwendung – ein Traum für unser Stromnetz und unsere Umwelt. Oder Züge, die auf einem Magnetfeld schweben und lautlos durch die Landschaft gleiten. All das ist die Vision, die die Forschung an Supraleitern antreibt.

Was sind Supraleiter? Ein kurzer Blick

Supraleiter sind Materialien, die bei extrem niedrigen Temperaturen (nahe dem absoluten Nullpunkt von -273.15 °C) zwei bemerkenswerte Eigenschaften aufweisen:

  1. Null elektrischer Widerstand: Strom fließt durch sie ohne jegliche Energieverluste.
  2. Meißner-Effekt: Sie stoßen Magnetfelder vollständig ab, was zum berühmten „Schweben“ von Magneten über dem Material führt.

Lange Zeit waren Supraleiter reine Tieftemperaturphänomene, beschränkt auf wenige Materialien und Kühlanwendungen mit flüssigem Helium. Doch dann, 1986, kam die Revolution: Die Entdeckung der Hochtemperatur-Supraleiter (HTSL). Diese komplexen Keramiken, oft auf Kupferoxid-Basis (sogenannte Kuprate), konnten bei deutlich höheren Temperaturen (immer noch kalt, aber erreichbar mit flüssigem Stickstoff bei -196 °C) supraleitend werden. Dies öffnete die Tür zu potenziellen kommerziellen Anwendungen, auch wenn der breite Einsatz noch vor großen technologischen Hürden steht.


Das Herz der Kuprate: Die CuO2-Ebene

Im Zentrum fast aller Hochtemperatur-Supraleiter auf Kuprat-Basis steht eine atomar dünne Schicht: die Kupferoxid-Ebene (CuO2-Ebene).

Diese Ebene besteht aus Kupferatomen, die in einem quadratischen Gitter angeordnet sind, wobei jedes Kupferatom von vier Sauerstoffatomen umgeben ist und ein Kreuz bildet. Man kann sich diese Ebenen als die „Leitungsbahnen“ vorstellen, in denen die faszinierenden elektronischen Phänomene der Supraleitung stattfinden. Die komplexen Theorien zur Hochtemperatur-Supraleitung drehen sich oft darum, wie Elektronen und ihre Spins in dieser speziellen Umgebung miteinander wechselwirken.


Die Sprache der Elektronen: Hubbard-Parameter aus „First Principles“

Um das Verhalten von Elektronen in diesen komplexen Materialien zu beschreiben, verwenden Physiker oft vereinfachte, aber leistungsstarke Modelle. Das Hubbard-Modell ist ein solches fundamentales Modell, das die wichtigsten Interaktionen zwischen Elektronen auf einem Gitter erfasst. Es basiert auf wenigen Schlüsselparametern, die die elektronische Struktur und das magnetische Verhalten bestimmen:

  • Ud​ (On-site Coulomb-Abstoßung): Dies beschreibt die Energie, die nötig ist, um zwei Elektronen auf demselben Atomorbital (hier: Kupfer 3d-Orbital) zu platzieren. Es ist eine Maß für die Abstoßung zwischen Elektronen. Ein hoher Ud führt dazu, dass Elektronen es vorziehen, auf verschiedenen Atomen zu sitzen.
  • tpd (Hopping-Integral): Dies beschreibt die Wahrscheinlichkeit, dass ein Elektron von einem Sauerstoff-2p-Orbital auf ein benachbartes Kupfer-3d-Orbital springt (oder „hüpft“). Es ist ein Maß für die Delokalisierung der Elektronen und damit für die kinetische Energie.
  • Δ (Ladungstransferenergie): Dies ist die Energie, die benötigt wird, um ein Elektron von einem Sauerstoff-2p-Orbital auf ein Kupfer-3d-Orbital zu transferieren. Man kann es sich als die Energiedifferenz zwischen den lokalen p- und d-Niveaus vorstellen.

Unsere Berechnung am CuO2-Cluster

Traditionell wurden diese Parameter oft aus experimentellen Daten abgeleitet. Doch moderne Quantenchemie-Methoden aus „First Principles“ (d.h., nur basierend auf den fundamentalen physikalischen Konstanten und den Atomkernen, ohne experimentelle Anpassung) erlauben es uns, diese Parameter direkt zu berechnen.

Um die Parameter für die CuO2-Ebene zu erhalten, haben wir nicht die gesamte unendliche Schicht simuliert (was extrem rechenintensiv wäre), sondern ein repräsentatives kleines Fragment: ein CuO2-Cluster mit einem 90-Grad-Winkel. Dieses Cluster fängt die lokale Umgebung eines Kupferatoms in der Ebene ein.

Wir haben die PySCF-Bibliothek für unsere quantenchemischen Berechnungen verwendet.

Berechnung von Ud​ (On-site Coulomb-Abstoßung): Für Ud​ betrachten wir ein isoliertes Kupferatom in verschiedenen Ladungszuständen (Cu+,Cu2+,Cu3+). Die Energiedifferenzen zwischen diesen Zuständen liefern uns einen Wert für die Abstoßung. Resultat: Ud​≈17.53 eV

Berechnung von tpd (Hopping-Integral): Wir platzieren das Kupferatom in der Mitte des Clusters und die beiden Sauerstoffatome im Abstand von 1.9 Ångström. Die „Hopping“-Stärke tpd​ wird direkt aus dem Core-Hamiltonian der berechneten elektronischen Struktur extrahiert, der die Wechselwirkung zwischen dem Kupfer 3dx2-y2-Orbital und den Sauerstoff 2p-Orbitalen beschreibt. Resultat: tpd≈3.18 eV (Betrag)

Berechnung von Δ (Ladungstransferenergie): Δ wurde als Energiedifferenz zwischen dem höchsten besetzten Molekülorbital (HOMO) mit dominantem Sauerstoff-2p-Charakter und dem niedrigsten unbesetzten Molekülorbital (LUMO) mit dominantem Kupfer-3d-Charakter im CuO2-Cluster abgeschätzt. Resultat: Δ≈13.65 eV

Es ist wichtig zu beachten, dass diese „nackten“ Parameter aus kleinen Clustern und unter Vernachlässigung von Abschirmungseffekten (die im Festkörper auftreten würden) in der Regel von den „effektiven“ Parametern abweichen, die in detaillierten Festkörpermodellen verwendet werden.


Der Sprung zu „J“: Magnetismus in Kupraten

Warum sind diese Parameter so wichtig? Sie fließen in das t-J-Modell ein, ein fundamentales Modell zur Beschreibung der elektronischen Struktur und des Magnetismus in Kupraten.

Das t-J-Modell ist eine Vereinfachung des Hubbard-Modells, das besonders gut die starken Elektron-Elektron-Wechselwirkungen in Kupraten erfasst. Eine seiner Schlüsselgrößen ist die magnetische Austauschkopplung J. Dieses J beschreibt die Stärke der Wechselwirkung zwischen den Spins benachbarter Kupferatome. In Kupraten ist J typischerweise antiferromagnetisch, was bedeutet, dass benachbarte Spins sich bevorzugt antiparallel ausrichten.

Die Beziehung zwischen unseren berechneten Parametern und J wird durch die Anderson-Superexchange-Formel gegeben:

J≈4tpd4​​ / Δ2Ud

Mit unseren Werten:

  • tpd≈3.18 eV
  • Ud​≈17.53 eV
  • Δ≈13.65 eV

erhalten wir:

J≈4(3.18 eV)4 /((13.65 eV)2⋅(17.53 eV))≈0.1251 eV

Das ist ein außergewöhnlich gutes Ergebnis! Experimentelle Werte für J in Hochtemperatur-Supraleitern liegen typischerweise im Bereich von 0.1 bis 0.15 eV. Unser berechneter Wert passt erstaunlich gut in dieses Fenster, auch wenn wir nur einen kleinen Cluster und „nackte“ Parameter verwendet haben. Dies deutet darauf hin, dass die grundlegenden physikalischen Wechselwirkungen innerhalb dieser CuO2-Cluster gut erfasst werden.


Das t-J-Modell: Ein Kandidat für die HTC-Theorie

Das t-J-Modell ist besonders attraktiv für die Beschreibung von Hochtemperatur-Supraleitern, weil es die zwei wichtigsten Aspekte dieser Materialien vereint:

  1. Starke Korrelationen: Es berücksichtigt die starke Coulomb-Abstoßung (Ud​), die verhindert, dass Elektronen dasselbe Atom besetzen, und führt zu sogenannten stark korrelierten Systemen.
  2. Spin-Dynamik: Es enthält explizit die magnetische Austauschkopplung (J), die für die antiferromagnetische Ordnung der Spins der Kupferatome verantwortlich ist.

Der Glaube ist, dass die Supraleitung in Kupraten durch die Dotierung dieser antiferromagnetisch geordneten Ebenen entsteht. Wenn man Ladungsträger (Löcher) in diese Ebene einbringt, stören sie die magnetische Ordnung. Die Wechselwirkung dieser Ladungsträger mit den Spin-Fluktuationen, die durch J vermittelt werden, könnte der „Klebstoff“ sein, der Elektronen zu Cooper-Paaren bindet und sie dann supraleitend macht. Das t-J-Modell bietet einen Rahmen, um diese komplexen Wechselwirkungen zu studieren.

Obwohl das t-J-Modell eine vereinfachte Darstellung ist, hat es sich als äußerst nützlich erwiesen, um viele Eigenschaften von HTSL zu erklären, einschließlich ihrer einzigartigen Normalzustands-Eigenschaften (oberhalb der Sprungtemperatur) und der Supraleitung selbst. Unsere „First-Principles“-Berechnungen leisten einen Beitrag dazu, die Parameter dieses Modells direkt aus der fundamentalen Quantenmechanik abzuleiten und so die Brücke zwischen ab-initio-Berechnungen und phänomenologischen Modellen zu schlagen.

Source Code

Eine grobe Abschätzung für Ud:

import numpy
from pyscf import gto, scf

# Define a single Copper atom at the origin
atom_coords = [['Cu', (0, 0, 0)]]
basis_set = 'def2-svp' # Use the same basis set as for Cu-O

# --- 1. Berechnung für Cu+ (d10 Konfiguration) ---
# Copper has 29 electrons. Cu+ has 28 electrons.
# This is a closed-shell system (all electrons paired).
mol_cu_plus = gto.M(
    atom = atom_coords,
    basis = basis_set,
    charge = 1, # +1 charge
    spin = 0,   # Closed shell (28 electrons)
    verbose = 0, # Less verbose for cleaner output
)
mf_cu_plus = scf.RHF(mol_cu_plus).run()
E_cu_plus = mf_cu_plus.e_tot
print(f"Total energy of Cu+ (d10): {E_cu_plus:.8f} Hartree")

# --- 2. Berechnung für Cu2+ (d9 Konfiguration) ---
# Copper has 29 electrons. Cu2+ has 27 electrons.
# This is an open-shell system (odd number of electrons, S=1/2).
# We use UHF (Unrestricted Hartree-Fock) for open shells.
mol_cu_2plus = gto.M(
    atom = atom_coords,
    basis = basis_set,
    charge = 2, # +2 charge
    spin = 1,   # 2S = 1 for S=1/2 (one unpaired electron)
    verbose = 0,
)
mf_cu_2plus = scf.UHF(mol_cu_2plus).run() # Use UHF for open shells
E_cu_2plus = mf_cu_2plus.e_tot
print(f"Total energy of Cu2+ (d9): {E_cu_2plus:.8f} Hartree")

# --- 3. Berechnung für Cu3+ (d8 Konfiguration) ---
# Copper has 29 electrons. Cu3+ has 26 electrons.
# This is an open-shell system (even number of electrons, but can be singlet S=0 or triplet S=1).
# For d8, the lowest energy state is often a triplet (spin=2). Let's calculate for triplet first.
mol_cu_3plus_triplet = gto.M(
    atom = atom_coords,
    basis = basis_set,
    charge = 3, # +3 charge
    spin = 2,   # 2S = 2 for S=1 (two unpaired electrons)
    verbose = 0,
)
mf_cu_3plus_triplet = scf.UHF(mol_cu_3plus_triplet).run()
E_cu_3plus_triplet = mf_cu_3plus_triplet.e_tot
print(f"Total energy of Cu3+ (d8, Triplet): {E_cu_3plus_triplet:.8f} Hartree")

# Optionally, check singlet for Cu3+ as well (often higher in energy)
# mol_cu_3plus_singlet = gto.M(
#     atom = atom_coords,
#     basis = basis_set,
#     charge = 3, # +3 charge
#     spin = 0,   # 2S = 0 for S=0 (all paired)
#     verbose = 0,
# )
# mf_cu_3plus_singlet = scf.UHF(mol_cu_3plus_singlet).run() # UHF can optimize to singlet if it's lowest
# E_cu_3plus_singlet = mf_cu_3plus_singlet.e_tot
# print(f"Total energy of Cu3+ (d8, Singlet): {E_cu_3plus_singlet:.8f} Hartree")

# Choose the lowest energy for Cu3+
E_cu_3plus = E_cu_3plus_triplet # Assuming triplet is lower for d8


# --- 4. Berechnung von U_d ---
# U_d = E(Cu3+) + E(Cu+) - 2 * E(Cu2+)
U_d_hartree = E_cu_3plus + E_cu_plus - 2 * E_cu_2plus
U_d_eV = U_d_hartree * 27.2114

print(f"\nCalculated U_d for isolated Cu atom: {U_d_hartree:.4f} Hartree")
print(f"                                   : {U_d_eV:.4f} eV")

# --- Vergleich mit Literaturwerten ---
print("\n--- Comparison with literature values for Hubbard U ---")
print("Literature values for U_d in solid cuprates are typically around 8-10 eV (screened).")
print("The calculated value for an isolated atom is expected to be higher due to lack of screening.")

Eine grobe Abschätzung für tpd und Δ:

import numpy
from pyscf import gto, scf

# --- 1. System Geometrie und Basissätze definieren ---
cu_o_distance = 1.9 # Ångström, typical Cu-O distance

mol = gto.M(
    atom = [
        ['Cu', (0, 0, 0)],           # Atom 0: Central Copper
        ['O',  (cu_o_distance, 0, 0)], # Atom 1: Oxygen on X-axis
        ['O',  (0, cu_o_distance, 0)], # Atom 2: Oxygen on Y-axis
    ],
    basis = {
        'Cu': 'def2-svp',
        'O':  'def2-svp',
    },
    charge = 0, # Neutral system
    spin = 1,   # Doublet state (one unpaired electron on Cu 3d)
    verbose = 4, # Keep verbose for detailed output
)

# --- 2. SCF Berechnung (UHF für Open-Shell-Systeme) ---
mf = scf.UHF(mol).run()

# --- 3. Extrahiere Ein-Elektronen-Integrale (Core Hamiltonian) ---
hcore = mf.get_hcore()

print("\n--- BEGIN FULL AO LABELS LIST (copied from mol.ao_labels()) ---")
ao_labels = mol.ao_labels()
for i, label in enumerate(ao_labels):
    print(f"Index {i}: '{label}'") # Print with quotes to make hidden spaces visible
print("--- END FULL AO LABELS LIST ---\n")

# --- 4. Identifiziere relevante Orbitale für t_pd ---
# Dies muss sehr genau auf die AO-Label-Strings von PySCF abgestimmt sein.
# Überprüfen Sie Ihre Ausgabe von mol.ao_labels() erneut!
cu_d_x2y2_idx = -1
o1_px_idx = -1 # Oxygen 1 (on X-axis), 2px orbital
o2_py_idx = -1 # Oxygen 2 (on Y-axis), 2py orbital

for i, label in enumerate(ao_labels):
    # Cu 3d(x^2-y^2) - Atom 0
    if (label == '0 Cu 3dx2-y2 ') or (label.strip() == '0 Cu 3dx2-y2'):
        cu_d_x2y2_idx = i
    
    # O1 (Atom 1) 2px (sigma-bonding along X-axis)
    if (label == '1 O 2px   ') or (label.strip() == '1 O 2px'):
        o1_px_idx = i
        
    # O2 (Atom 2) 2py (sigma-bonding along Y-axis)
    if (label == '2 O 2py   ') or (label.strip() == '2 O 2py'):
        o2_py_idx = i

print(f"\nIdentified Indices for t_pd calculation:")
print(f"Central Cu 3d_x2-y2 index: {cu_d_x2y2_idx}")
print(f"O1 (on X-axis) 2p_x index: {o1_px_idx}")
print(f"O2 (on Y-axis) 2p_y index: {o2_py_idx}")


# --- 5. Berechne Hopping-Integrale t_pd ---
if cu_d_x2y2_idx == -1 or o1_px_idx == -1 or o2_py_idx == -1:
    print("\nError: Could not find all necessary orbital indices for t_pd. Check basis set or labels again.")
    t_pd_estimate_eV_avg = None # Initialize as None if not found
else:
    t_pd_hartree_O1 = hcore[cu_d_x2y2_idx, o1_px_idx]
    t_pd_eV_O1 = t_pd_hartree_O1 * 27.2114

    t_pd_hartree_O2 = hcore[cu_d_x2y2_idx, o2_py_idx]
    t_pd_eV_O2 = t_pd_hartree_O2 * 27.2114

    # --- FIX FOR t_pd AVERAGE ---
    # Take the absolute value of one of the t_pd values, as they are symmetric
    # Using the average of absolute values is also fine, but just one is simpler
    t_pd_estimate_eV_avg = abs(t_pd_eV_O1) # OR abs(t_pd_eV_O2), they should be the same magnitude

    print("\n--- Estimated Cu-O Hopping Integral (t_pd) in CuO2 cluster ---")
    print(f"t_pd (Cu-O1_px): {t_pd_hartree_O1:.4f} Hartree / {t_pd_eV_O1:.4f} eV")
    print(f"t_pd (Cu-O2_py): {t_pd_hartree_O2:.4f} Hartree / {t_pd_eV_O2:.4f} eV")
    print(f"Effective t_pd used for J (absolute value): {t_pd_estimate_eV_avg:.4f} eV")


# --- 6. Molecular Orbital Analysis for Delta Estimation ---
print("\n--- Molecular Orbital Analysis for Delta Estimation ---")

# MO-Energien und Koeffizienten für Alpha-Spin
mo_energies_alpha = mf.mo_energy[0]
mo_coeffs_alpha = mf.mo_coeff[0]

# Anzahl der besetzten Alpha-Orbitale
nocc_alpha = mol.nelectron // 2 + mol.spin # For a doublet, this is (N_total + 1) / 2
print(f"Number of alpha occupied orbitals: {nocc_alpha}")
print(f"Total number of alpha MOs: {len(mo_energies_alpha)}")

# Identifiziere die Indizes für *alle* Cu 3d und O 2p Atomorbitale
cu_3d_ao_indices = []
o_2p_ao_indices = []

for i, label in enumerate(ao_labels):
    # Cu 3d orbitals (all 5 of them) - Atom 0
    if ('0 Cu 3dxy ' == label) or ('0 Cu 3dxy' == label.strip()): cu_3d_ao_indices.append(i)
    elif ('0 Cu 3dyz ' == label) or ('0 Cu 3dyz' == label.strip()): cu_3d_ao_indices.append(i)
    elif ('0 Cu 3dz^2 ' == label) or ('0 Cu 3dz^2' == label.strip()): cu_3d_ao_indices.append(i)
    elif ('0 Cu 3dxz ' == label) or ('0 Cu 3dxz' == label.strip()): cu_3d_ao_indices.append(i)
    elif ('0 Cu 3dx2-y2 ' == label) or ('0 Cu 3dx2-y2' == label.strip()): cu_3d_ao_indices.append(i)

    # O1 2p orbitals (all 3 of them) - Atom 1
    if ('1 O 2px   ' == label) or ('1 O 2px' == label.strip()): o_2p_ao_indices.append(i)
    elif ('1 O 2py   ' == label) or ('1 O 2py' == label.strip()): o_2p_ao_indices.append(i)
    elif ('1 O 2pz   ' == label) or ('1 O 2pz' == label.strip()): o_2p_ao_indices.append(i)

    # O2 2p orbitals (all 3 of them) - Atom 2
    if ('2 O 2px   ' == label) or ('2 O 2px' == label.strip()): o_2p_ao_indices.append(i)
    elif ('2 O 2py   ' == label) or ('2 O 2py' == label.strip()): o_2p_ao_indices.append(i)
    elif ('2 O 2pz   ' == label) or ('2 O 2pz' == label.strip()): o_2p_ao_indices.append(i)


print(f"Indices of Cu 3d AOs for MO character analysis: {cu_3d_ao_indices}")
print(f"Indices of O 2p AOs for MO character analysis: {o_2p_ao_indices}")


# --- Identifiziere O 2p-ähnliches HOMO und Cu 3d-ähnliches LUMO ---
o_p_homo_energy = -float('inf')
o_p_homo_idx = -1
if o_2p_ao_indices:
    for i in range(nocc_alpha - 1, -1, -1): # Iterate from HOMO (nocc_alpha-1) downwards
        mo_contribution_o_p = 0.0
        for ao_idx in o_2p_ao_indices:
            mo_contribution_o_p += mo_coeffs_alpha[ao_idx, i]**2

        if mo_contribution_o_p > 0.30: # Threshold for O 2p character
            if mo_energies_alpha[i] > o_p_homo_energy:
                o_p_homo_energy = mo_energies_alpha[i]
                o_p_homo_idx = i
else:
    print("Warning: No O 2p AOs found for character analysis.")


# Find the lowest energy MO with significant Cu 3d character among unoccupied orbitals
cu_d_lumo_energy = float('inf')
cu_d_lumo_idx = -1
if cu_3d_ao_indices:
    print("\n--- Unoccupied MOs with Cu 3d Character (for Delta diagnostics) ---")
    for i in range(nocc_alpha, len(mo_energies_alpha)): # Iterate from LUMO (nocc_alpha) upwards
        mo_contribution_cu_d = 0.0
        for ao_idx in cu_3d_ao_indices:
            mo_contribution_cu_d += mo_coeffs_alpha[ao_idx, i]**2

        print(f"MO {i}: Energy: {mo_energies_alpha[i]:.4f} Hartree, Cu 3d Char: {mo_contribution_cu_d*100:.1f}%")

        if mo_contribution_cu_d > 0.10: # Threshold for Cu 3d character
            if mo_energies_alpha[i] < cu_d_lumo_energy:
                cu_d_lumo_energy = mo_energies_alpha[i]
                cu_d_lumo_idx = i
            # --- FIX: REMOVE THE 'break' STATEMENT HERE ---
            # break # <--- REMOVE THIS LINE!
else:
    print("Warning: No Cu 3d AOs found for character analysis.")


print(f"\nEstimated O 2p-like HOMO (alpha): MO {o_p_homo_idx}, Energy: {o_p_homo_energy:.4f} Hartree")
print(f"Estimated Cu 3d-like LUMO (alpha): MO {cu_d_lumo_idx}, Energy: {cu_d_lumo_energy:.4f} Hartree")

delta_eV = None
if o_p_homo_idx != -1 and cu_d_lumo_idx != -1:
    delta_hartree = cu_d_lumo_energy - o_p_homo_energy
    delta_eV = delta_hartree * 27.2114
    print(f"Rough estimated Delta (charge transfer energy): {delta_hartree:.4f} Hartree")
    print(f"                                             : {delta_eV:.4f} eV")
else:
    print("Could not find suitable O 2p-like HOMO or Cu 3d-like LUMO for Delta estimation.")
    print("Adjust the character thresholds (e.g., 0.30) or print more MO data manually for inspection.")


# --- 7. Abschätzung von J mit den berechneten "Bare"-Werten ---
# t_pd_estimate_eV_avg from this calculation
u_d_eV = 17.5255 # Your previously calculated U_d for isolated Cu atom

if 't_pd_estimate_eV_avg' in locals() and t_pd_estimate_eV_avg is not None and delta_eV is not None:
    if abs(delta_eV) < 1e-6: # Prevent division by zero
        J_estimated_eV = float('inf')
        print("\nWarning: Estimated Delta is too close to zero, J calculation might be unstable.")
    else:
        J_estimated_eV = (4 * (t_pd_estimate_eV_avg**4)) / ((delta_eV**2) * u_d_eV)
        
    print(f"\n--- Estimated J with bare t_pd, bare U_d, and rough Delta ---")
    print(f"J = 4 * (t_pd^4) / (Delta^2 * U_d)")
    print(f"J approx 4 * ({t_pd_estimate_eV_avg:.4f}^4) / ({delta_eV:.4f}^2 * {u_d_eV:.4f})")
    print(f"J estimated: {J_estimated_eV:.4f} eV")
    print("\n--- Comparison with literature values ---")
    print("Experimental J in solid cuprates: ~0.1 - 0.15 eV.")
    print("This calculated J will likely be much higher or lower due to using unscreened parameters from a small cluster.")
else:
    print("\nCannot estimate J because t_pd_average or Delta could not be determined.")