Skip to content

Commit eebf8d8

Browse files
committed
additional notes
1 parent 18ec7a1 commit eebf8d8

File tree

3 files changed

+219
-0
lines changed

3 files changed

+219
-0
lines changed

content/notes/2026-02-PME.md

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
---
2+
title: Smooth Particle Mesh Ewald (PME)
3+
date: 2025-02-15
4+
math: true
5+
6+
image:
7+
filename: ""
8+
caption: "Math / formula (optional caption)"
9+
focal_point: Smart
10+
preview_only: false
11+
---
12+
13+
# Smooth Particle Mesh Ewald (PME)
14+
15+
In molecular simulations, calculating electrostatic interactions is computationally expensive because the Coulomb potential ($1/r$) is long-ranged. A "brute-force" calculation for $N$ particles scales as $O(N^2)$. **Smooth Particle Mesh Ewald (PME)** is an algorithm that reduces this complexity to $O(N \log N)$.
16+
17+
## 1. The Core Problem
18+
When using periodic boundary conditions, we must sum the interactions of a particle with every other particle in an infinite periodic grid. This sum is only conditionally convergent and extremely slow to compute directly.
19+
20+
The Ewald technique splits the interaction into two parts:
21+
1. **Short-range (Real Space):** Interactions that decay quickly.
22+
2. **Long-range (Reciprocal Space):** A slowly varying part handled in Fourier space.
23+
24+
---
25+
26+
## 2. Mathematical Decomposition
27+
28+
The total electrostatic energy $V_{total}$ is expressed as:
29+
30+
$$V_{total} = V_{real} + V_{recip} + V_{self}$$
31+
32+
### A. Real-Space Term (Short Range)
33+
This term handles nearby particles. It uses the **complementary error function** ($\text{erfc}$) to screen the charges:
34+
35+
$$V_{real} = \frac{1}{2} \sum_{i,j}^N \sum_{n \in \mathbb{Z}^3}^* \frac{q_i q_j \text{erfc}(\alpha |r_{ij} + nL|)}{|r_{ij} + nL|}$$
36+
37+
* **$\alpha$**: The Ewald splitting parameter (controls the rate of decay).
38+
* **$nL$**: Periodic images in a box of length $L$.
39+
40+
41+
42+
### B. Reciprocal-Space Term (Long Range)
43+
The "smooth" part of the potential is periodic and solved using a Fourier Series:
44+
45+
$$V_{recip} = \frac{1}{2\pi \Omega} \sum_{k \neq 0} \frac{4\pi}{k^2} e^{-k^2 / 4\alpha^2} |S(k)|^2$$
46+
47+
* **$\Omega$**: The volume of the unit cell.
48+
* **$S(k)$**: The **structure factor**, $\sum q_j e^{ik \cdot r_j}$.
49+
50+
> **The PME Innovation:** Instead of calculating $S(k)$ directly, PME interpolates charges onto a 3D grid using **B-splines** and solves the sum using the **Fast Fourier Transform (FFT)**.
51+
52+
### C. Self-Correction Term
53+
Since the reciprocal sum includes the interaction of a particle with its own neutralizing background, we subtract the self-energy:
54+
55+
$$V_{self} = -\frac{\alpha}{\sqrt{\pi}} \sum_{i=1}^N q_i^2$$
56+
57+
---
58+
59+
## 3. The PME Workflow
60+
61+
62+
63+
1. **Charge Assignment:** Map point charges $q_i$ to a regular 3D mesh using cardinal B-splines.
64+
2. **FFT:** Transform the grid-based charge distribution into reciprocal space.
65+
3. **Green's Function:** Multiply by the reciprocal potential (the $1/k^2$ term).
66+
4. **Inverse FFT:** Transform back to real space to get the potential on the grid.
67+
5. **Force Interpolation:** Calculate the gradient of the potential at each particle's actual position to determine the forces.
68+
69+
---
70+
71+
## 4. Key Parameters in MD Engines
72+
73+
| Parameter | Recommended Value / Impact |
74+
| :--- | :--- |
75+
| **PME Order** | Usually **4** (Cubic Splines). Balance between accuracy and speed. |
76+
| **Grid Spacing** | Typically **0.10 to 0.15 nm**. Affects FFT resolution. |
77+
| **$\alpha$ (Tolerance)** | Often derived from the real-space cutoff (e.g., $10^{-5}$). |
78+
79+
80+
81+
82+
83+
---
84+
85+
## 5. Foundational & Modern References
86+
87+
### The Smooth Particle Mesh Ewald (SPME)
88+
* **Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., & Pedersen, L. G. (1995).** "A smooth particle mesh Ewald method." *The Journal of Chemical Physics*, 103(19), 8577-8593.
89+
90+
### The Original Particle Mesh Ewald (PME)
91+
* **Darden, T., York, D., & Pedersen, L. (1993).** "Particle mesh Ewald: An $N \cdot \log(N)$ method for Ewald sums in large systems." *The Journal of Chemical Physics*, 98(12), 10089-10092.
92+
93+
### The Original Ewald Method
94+
* **Ewald, P. P. (1921).** "Die Berechnung optischer und elektrostatischer Gitterpotentiale." *Annalen der Physik*, 369(3), 253-287.
95+
96+
97+
### Modern Implementation: DIMOS (PyTorch)
98+
* **Christiansen, H., et al. (2025).** "Fast, Modular, and Differentiable Framework for Machine Learning-Enhanced Molecular Simulations." *The Journal of Chemical Physics*.
99+
* [Article DOI: 10.1063/5.0277356](https://doi.org/10.1063/5.0277356)
100+
* [arXiv: 2503.20541](https://arxiv.org/abs/2503.20541)
101+
* [GitHub Code: nec-research/DIMOS](https://github.com/nec-research/DIMOS)
102+
> **Note:** This framework provides a high-performance, differentiable environment for integrating ML potentials directly into MD/MC workflows using PyTorch.
103+
104+
105+
106+
---
107+
*Note: PME is essential for maintaining stability in simulations of charged systems like DNA, proteins, and lipid bilayers.*
108+
109+
---
110+
*Note: This note was prepared as a technical summary for molecular dynamics enthusiasts.*
Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
---
2+
title: Spherical Harmonics
3+
date: 2025-02-15
4+
math: true
5+
image:
6+
filename: ""
7+
caption: ""
8+
focal_point: Smart
9+
preview_only: false
10+
---
11+
12+
# Spherical Harmonics and e3nn
13+
14+
**Spherical Harmonics** ($Y_{\ell}^m$) are a set of orthogonal functions defined on the surface of a sphere. They are the angular portion of a set of solutions to Laplace's equation and serve as the Fourier basis for functions defined on a sphere.
15+
16+
## 1. Mathematical Overview
17+
18+
In spherical coordinates $(\theta, \phi)$, the real spherical harmonics are defined as:
19+
$$Y_{\ell}^m(\theta, \phi) = N_{\ell}^m P_{\ell}^m(\cos \theta) \text{trig}(m\phi)$$
20+
21+
Where:
22+
- $\ell \ge 0$ is the **degree** (angular momentum).
23+
- $-\ell \le m \le \ell$ is the **order**.
24+
- $P_{\ell}^m$ are the associated Legendre polynomials.
25+
- $N_{\ell}^m$ is a normalization constant.
26+
27+
In modern geometric machine learning, we often represent these in **Cartesian coordinates** $(x, y, z)$, where they become homogeneous polynomials of degree $\ell$.
28+
29+
---
30+
31+
## 2. Using the `e3nn` Library
32+
33+
The [e3nn](https://github.com/e3nn/e3nn) library is designed for **Euclidean Neural Networks**. It uses spherical harmonics as the primary way to lift geometric data (like atomic positions) into a representation that is equivariant to rotations ($SO(3)$) and inversion ($O(3)$).
34+
35+
### Key Concepts in e3nn:
36+
* **Irreps (Irreducible Representations):** e3nn tracks how data transforms. A spherical harmonic of degree $\ell$ belongs to the irrep $\ell$.
37+
* **Parity:** In e3nn, spherical harmonics have a definite parity $p = (-1)^\ell$.
38+
* **Coordinate Convention:** e3nn typically expects input vectors in $(x, y, z)$ but internally uses a $y, z, x$ convention for alignment with standard real spherical harmonics.
39+
40+
---
41+
42+
## 3. Python Implementation
43+
44+
To use this code, you will need to install the library:
45+
`pip install e3nn torch`
46+
47+
### Code Example: Generating and Visualizing
48+
The following script calculates the spherical harmonic coefficients for a given set of vectors.
49+
50+
```python
51+
import torch
52+
from e3nn import o3
53+
import numpy as np
54+
import matplotlib.pyplot as plt
55+
from matplotlib import cm
56+
57+
def plot_all_harmonics(l_max):
58+
fig = plt.figure(figsize=(12, 8))
59+
60+
# Grid coordinates for sampling the sphere
61+
phi, theta = np.mgrid[0:2*np.pi:70j, 0:np.pi:70j]
62+
x = np.sin(theta) * np.cos(phi)
63+
y = np.sin(theta) * np.sin(phi)
64+
z = np.cos(theta)
65+
vectors = torch.tensor(np.stack([x.flatten(), y.flatten(), z.flatten()], axis=1), dtype=torch.float32)
66+
67+
plot_idx = 1
68+
for l in range(l_max + 1):
69+
# Get irreps and compute SH for this specific L
70+
irreps = o3.Irreps.spherical_harmonics(l)
71+
sh_values = o3.spherical_harmonics(irreps, vectors, normalize=True)
72+
73+
# There are 2L + 1 components for each L
74+
num_m = 2 * l + 1
75+
76+
for m_idx in range(num_m):
77+
# Create subplot: rows = l_max+1, cols = 2*l_max+1
78+
# We center the plots by offsetting the start index
79+
ax_idx = l * (2 * l_max + 1) + (l_max - l) + m_idx + 1
80+
ax = fig.add_subplot(l_max + 1, 2 * l_max + 1, ax_idx, projection='3d')
81+
82+
# Extract the specific m component
83+
r_vals = sh_values[:, m_idx].reshape(70, 70).numpy()
84+
85+
# Radial distance is absolute value, color is the sign
86+
R = np.abs(r_vals)
87+
X_p, Y_p, Z_p = R * x, R * y, R * z
88+
89+
# Normalize colors so 0 is white/middle
90+
v_max = np.max(np.abs(r_vals))
91+
norm = plt.Normalize(-v_max, v_max)
92+
colors = cm.RdBu(norm(r_vals))
93+
94+
ax.plot_surface(X_p, Y_p, Z_p, facecolors=colors, antialiased=True, shade=True)
95+
96+
ax.set_title(f"L={l}, m_idx={m_idx}", fontsize=8)
97+
ax.axis('off')
98+
99+
plt.tight_layout()
100+
plt.savefig('spherical_harmonics_l3.png', dpi=300, transparent=True)
101+
plt.show()
102+
103+
plot_all_harmonics(l_max=3)
104+
105+
```
106+
107+
108+
![Spherical Harmonics plot](spherical_harmonics_l3.png)
109+
560 KB
Loading

0 commit comments

Comments
 (0)