.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/advanced/06_ewald_electrostatics.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_advanced_06_ewald_electrostatics.py: Long-Range Electrostatics with Ewald Summation =============================================== Lennard-Jones and other short-range potentials truncate interactions at a cutoff radius, which works well for neutral systems with rapidly decaying forces. **Ionic systems** — salts, molten metals, charged proteins — require long-range Coulomb interactions that decay as 1/r and cannot be safely truncated without severe artefacts. The **Ewald summation** method handles this by splitting the Coulomb sum into: * A **real-space** term — short-range, erfc-damped, evaluated over a neighbor list exactly as for LJ. * A **reciprocal-space** term — smooth structure-factor sum in Fourier space that captures the long-range tail. Together they reproduce the exact infinite-lattice Coulomb energy for a periodic system. This example: * Builds a minimal NaCl rock-salt crystal (2×2×2 supercell, 64 atoms) with correct formal charges +1 e (Na) and −1 e (Cl). * Computes energy and forces with :class:`~nvalchemi.models.ewald.EwaldModelWrapper`. * Compares the total Ewald energy to the analytical Madelung energy. .. note:: For molecular-dynamics simulations of ionic systems you need a short-range repulsion term alongside the Ewald electrostatics. Pure Coulomb has no equilibrium — without repulsion ions accelerate toward each other without bound. See **Example 07** (Additive Model Composition) for the complete ``LennardJones + Ewald`` setup used in production MD runs. Key concepts demonstrated -------------------------- * Constructing an :class:`~nvalchemi.data.AtomicData` with ``node_charges`` (shape ``[N, 1]``, elementary charge units). * Instantiating :class:`~nvalchemi.models.ewald.EwaldModelWrapper` with a real-space cutoff and auto-estimated Ewald parameters. * Computing energy and forces for a single snapshot via a direct ``model(batch)`` call. * Using :meth:`~nvalchemi.models.ewald.EwaldModelWrapper.invalidate_cache` when the simulation cell changes between evaluations. .. GENERATED FROM PYTHON SOURCE LINES 61-70 .. code-block:: Python from __future__ import annotations import torch from nvalchemi.data import AtomicData, Batch from nvalchemi.dynamics.hooks import NeighborListHook from nvalchemi.models.ewald import EwaldModelWrapper .. GENERATED FROM PYTHON SOURCE LINES 71-78 Ewald model setup ------------------ :class:`~nvalchemi.models.ewald.EwaldModelWrapper` only needs a real-space cutoff. The Ewald splitting parameter α and reciprocal-space cutoff are estimated automatically from the ``accuracy`` target (default 1e-6) each time the simulation cell changes. For a fixed-cell run the cache is computed once and reused. .. GENERATED FROM PYTHON SOURCE LINES 78-83 .. code-block:: Python CUTOFF = 10.0 # Å — real-space cutoff model = EwaldModelWrapper(cutoff=CUTOFF, accuracy=1e-6, max_neighbors=256) model.model_config.compute_forces = True .. GENERATED FROM PYTHON SOURCE LINES 84-97 Building a NaCl rock-salt supercell ------------------------------------- The NaCl rock-salt lattice has two interpenetrating FCC sub-lattices. We build a 2×2×2 conventional cubic supercell (64 atoms). **Important**: Na and Cl positions are collected separately, then concatenated, so that ``atomic_numbers`` and ``node_charges`` align correctly with ``positions``. Interleaving the two species within the same image ordering would silently mis-assign charges. Lattice constant a₀ = 5.64 Å (experimental value). Formal charges: Na → +1 e, Cl → −1 e. The system is overall charge-neutral (32 Na + 32 Cl). .. GENERATED FROM PYTHON SOURCE LINES 97-168 .. code-block:: Python A0 = 5.64 # Å — NaCl lattice constant SUPERCELL = 2 # repeat units along each axis M_NA = 22.990 # amu M_CL = 35.453 # amu cell_size = SUPERCELL * A0 # 11.28 Å per side # Fractional coordinates of the two sub-lattices within one conventional cell. # Na occupies the FCC vertices; Cl is offset by (0.5, 0, 0) etc. _na_basis = torch.tensor( [[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]], dtype=torch.float32, ) # 4 Na per unit cell _cl_basis = torch.tensor( [[0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5], [0.5, 0.5, 0.5]], dtype=torch.float32, ) # 4 Cl per unit cell # Supercell image offsets (in unit-cell units). offsets = torch.tensor( [ [ix, iy, iz] for ix in range(SUPERCELL) for iy in range(SUPERCELL) for iz in range(SUPERCELL) ], dtype=torch.float32, ) # (8, 3) # Build Na positions: (n_images, n_basis, 3) → (N_Na, 3). # Dividing by SUPERCELL converts unit-cell fractional to supercell fractional; # multiplying by cell_size converts to Cartesian Å. na_positions = ((_na_basis.unsqueeze(0) + offsets.unsqueeze(1)) / SUPERCELL).reshape( -1, 3 ) * cell_size cl_positions = ((_cl_basis.unsqueeze(0) + offsets.unsqueeze(1)) / SUPERCELL).reshape( -1, 3 ) * cell_size positions = torch.cat([na_positions, cl_positions], dim=0) # (64, 3) n_na, n_cl = na_positions.shape[0], cl_positions.shape[0] # 32, 32 n_atoms = n_na + n_cl # 64 # Atomic numbers and formal charges — correctly aligned with positions. atomic_numbers = torch.cat( [ torch.full((n_na,), 11, dtype=torch.long), torch.full((n_cl,), 17, dtype=torch.long), ] ) charges = torch.cat([torch.ones(n_na, 1), -torch.ones(n_cl, 1)]) # (N, 1) cell = torch.eye(3).unsqueeze(0) * cell_size # (1, 3, 3) data = AtomicData( positions=positions, atomic_numbers=atomic_numbers, node_charges=charges, # (N, 1) — required shape for AtomicData forces=torch.zeros(n_atoms, 3), energies=torch.zeros(1, 1), cell=cell, pbc=torch.tensor([[True, True, True]]), ) batch = Batch.from_data_list([data]) print( f"System: {n_na} Na + {n_cl} Cl, box={cell_size:.2f} Å, " f"nearest Na–Cl distance={A0 / 2:.3f} Å" ) .. rst-class:: sphx-glr-script-out .. code-block:: none System: 32 Na + 32 Cl, box=11.28 Å, nearest Na–Cl distance=2.820 Å .. GENERATED FROM PYTHON SOURCE LINES 169-173 Building the neighbor list and evaluating energy + forces ---------------------------------------------------------- For a one-shot energy evaluation, build the neighbor list manually using :class:`~nvalchemi.dynamics.hooks.NeighborListHook` outside the dynamics loop. .. GENERATED FROM PYTHON SOURCE LINES 173-185 .. code-block:: Python nl_hook = NeighborListHook(model.model_card.neighbor_config) nl_hook(batch, None) # populates batch.neighbor_matrix / batch.num_neighbors result = model(batch) energy_eV = result["energies"].item() forces = result["forces"] # (N, 3) eV/Å print(f"\nEwald energy: {energy_eV:.4f} eV") print(f"Max force magnitude: {forces.norm(dim=-1).max().item():.4f} eV/Å") .. rst-class:: sphx-glr-script-out .. code-block:: none Warp DeprecationWarning: The symbol `warp.types.Any` will soon be removed from the public API. It can still be accessed from `warp._src.types.Any` but might be changed or removed without notice. Ewald energy: -285.5506 eV Max force magnitude: 0.0000 eV/Å .. GENERATED FROM PYTHON SOURCE LINES 186-198 Comparison with the analytical Madelung energy ------------------------------------------------ For an infinite NaCl rock-salt crystal the total Coulomb energy per formula unit is: E₀ = −A * k_e * q² / r_{Na–Cl} where A = 1.7476 is the Madelung constant for the rock-salt structure, k_e = 14.3996 eV·Å/e², q = 1 e, and r_{Na–Cl} = a₀/2. For a finite supercell the Ewald result converges to the thermodynamic limit to within the accuracy set in the model (1e-6 here). .. GENERATED FROM PYTHON SOURCE LINES 198-212 .. code-block:: Python MADELUNG_NaCL = 1.7476 K_E = 14.3996 # eV·Å/e² r_nn = A0 / 2 # nearest-neighbour Na–Cl distance n_formula_units = n_na # 32 (one per Na–Cl pair) E_madelung = -MADELUNG_NaCL * K_E * 1.0**2 / r_nn * n_formula_units print(f"\nAnalytical Madelung energy ({n_formula_units} f.u.): {E_madelung:.4f} eV") print(f"Ewald energy: {energy_eV:.4f} eV") print( f"Relative error: {abs(energy_eV - E_madelung) / abs(E_madelung):.2e}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Analytical Madelung energy (32 f.u.): -285.5573 eV Ewald energy: -285.5506 eV Relative error: 2.36e-05 .. GENERATED FROM PYTHON SOURCE LINES 213-218 Forces should vanish at the perfect-crystal geometry ----------------------------------------------------- At the ideal rock-salt lattice sites every atom sits at a potential minimum by symmetry, so the net force on each atom is zero. Finite cell and numerical accuracy leave a small residual. .. GENERATED FROM PYTHON SOURCE LINES 218-224 .. code-block:: Python mean_force = forces.norm(dim=-1).mean().item() print( f"\nMean force magnitude at ideal geometry: {mean_force:.4e} eV/Å (should be ≈0)" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Mean force magnitude at ideal geometry: 5.3939e-07 eV/Å (should be ≈0) .. GENERATED FROM PYTHON SOURCE LINES 225-230 Evaluating with a perturbed geometry -------------------------------------- Displace atoms by a small random amount and re-evaluate to show non-zero forces. Call :meth:`~nvalchemi.models.ewald.EwaldModelWrapper.invalidate_cache` if the cell changes between evaluations (here the cell is unchanged, so not required). .. GENERATED FROM PYTHON SOURCE LINES 230-253 .. code-block:: Python torch.manual_seed(7) perturbed_positions = positions + 0.05 * torch.randn_like(positions) data_pert = AtomicData( positions=perturbed_positions, atomic_numbers=atomic_numbers, node_charges=charges, forces=torch.zeros(n_atoms, 3), energies=torch.zeros(1, 1), cell=cell, pbc=torch.tensor([[True, True, True]]), ) batch_pert = Batch.from_data_list([data_pert]) nl_hook(batch_pert, None) result_pert = model(batch_pert) print("\nPerturbed geometry:") print(f" Ewald energy: {result_pert['energies'].item():.4f} eV") print( f" Max force magnitude: {result_pert['forces'].norm(dim=-1).max().item():.4f} eV/Å" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Perturbed geometry: Ewald energy: -285.5374 eV Max force magnitude: 0.4427 eV/Å .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.344 seconds) .. _sphx_glr_download_examples_advanced_06_ewald_electrostatics.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06_ewald_electrostatics.ipynb <06_ewald_electrostatics.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06_ewald_electrostatics.py <06_ewald_electrostatics.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06_ewald_electrostatics.zip <06_ewald_electrostatics.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_