Симуляція газу Леннарда-Джонса з нуля

1. Чому Леннард-Джонс

Реальні атоми та малі молекули слабко притягуються на помірній відстані (сили Ван-дер-Ваальса / дисперсійні сили Лондона) і люто відштовхуються, щойно їхні електронні хмари починають перекриватись. Потенціал Леннарда-Джонса (ЛД), запропонований Джоном Леннардом-Джонсом у 1924 році, охоплює обидва ефекти єдиною простою формулою, яку дешево обчислювати мільйони разів на кадр — саме тому він і лишається типовою тестовою системою для навчання, бенчмаркінгу та вивчення загальної фазової поведінки (газ, рідина, тверде тіло) без прив'язки до будь-якої конкретної реальної хімічної речовини.

2. Потенціал і сила ЛД

Форма 12-6 балансує крутий відштовхувальний член (перекриття електронів, показник 12) проти м'якшого притягувального члена (дисперсія, показник 6):

Мінімум U(r) лежить при r = 2^(1/6) σ ≈ 1.122σ — природній рівноважній відстані між двома ізольованими частинками. Нижче цієї відстані сила сильно відштовхувальна; вище — слабко притягувальна, спадаючи до нуля при r → ∞.

function ljForceMagnitude(r, epsilon = 1, sigma = 1) {
  const sr6  = Math.pow(sigma / r, 6);
  const sr12 = sr6 * sr6;
  return 24 * epsilon / r * (2 * sr12 - sr6); // додатне = відштовхування, від'ємне = притягання
}

3. Зведені (ЛД) одиниці

Коди МД майже ніколи не симулюють в одиницях СІ. Замість того вони встановлюють ε = σ = m = 1 (маса частинки) і виражають усе інше — час, температуру, тиск — як безрозмірну комбінацію цих величин. Це тримає числа поблизу одиниці (добре поводяться чисельно) і робить одну симуляцію застосовною до будь-якої реальної речовини простим масштабуванням.

4. Інтегрування швидкісним Верле

Швидкісний Верле — стандартний інтегратор МД: симплектичний (зберігає об'єм фазового простору), оборотний у часі, і потребує лише одного обчислення сили на крок, попри те, що є другого порядку точності — критично важливо, коли обчислення сили — найдорожча частина кожного кадру.

function velocityVerletStep(particles, forces, dt, computeForces) {
  // Крок 1: оновити позиції, використовуючи старі прискорення
  for (const p of particles) {
    p.x += p.vx * dt + 0.5 * p.ax * dt * dt;
    p.y += p.vy * dt + 0.5 * p.ay * dt * dt;
  }
  // Крок 2: перерахувати сили при нових позиціях
  const oldAx = particles.map(p => p.ax);
  const oldAy = particles.map(p => p.ay);
  computeForces(particles); // встановлює p.ax, p.ay для нової конфігурації

  // Крок 3: оновити швидкості середнім старого й нового прискорення
  particles.forEach((p, i) => {
    p.vx += 0.5 * (oldAx[i] + p.ax) * dt;
    p.vy += 0.5 * (oldAy[i] + p.ay) * dt;
  });
}

5. Періодичні межі та найближче зображення

Ящик з кількома сотнями частинок домінується поверхневими ефектами, якщо не прибрати поверхню повністю: періодичні граничні умови нескінченно замощують ящик симуляції в кожному напрямку, тож частинка, що виходить через правий край, знову входить зліва. Обчислення сил тоді використовують умову найближчого зображення — для кожної пари беруть ту періодичну копію другої частинки, що найближча.

function wrapPosition(p, boxSize) {
  p.x = ((p.x % boxSize) + boxSize) % boxSize;
  p.y = ((p.y % boxSize) + boxSize) % boxSize;
}

function minimumImageDelta(dx, dy, boxSize) {
  // зсунути дельту в (-boxSize/2, boxSize/2], щоб пари "бачили" найближчу періодичну копію
  if (dx >  boxSize / 2) dx -= boxSize;
  if (dx <= -boxSize / 2) dx += boxSize;
  if (dy >  boxSize / 2) dy -= boxSize;
  if (dy <= -boxSize / 2) dy += boxSize;
  return { dx, dy };
}

6. Радіус обрізання і списки сусідів

За межами кількох кратних σ притягання ЛД нехтовно мале. Обрізання потенціалу на радіусі обрізання — за традицією r_c = 2.5σ — перетворює цикл сил всіх пар O(N²) на такий, що перевіряє лише пари в межах r_c, а просторова сітка (список комірок) робить цю перевірку в середньому O(N) замість O(N²).

function buildCellList(particles, boxSize, cutoff) {
  const nCells = Math.max(1, Math.floor(boxSize / cutoff));
  const cellSize = boxSize / nCells;
  const cells = Array.from({ length: nCells * nCells }, () => []);
  particles.forEach((p, i) => {
    const cx = Math.floor(p.x / cellSize) % nCells;
    const cy = Math.floor(p.y / cellSize) % nCells;
    cells[cy * nCells + cx].push(i);
  });
  return { cells, nCells, cellSize }; // перевіряти лише блок 3×3 сусідніх комірок для кожної частинки
}

7. Термостатування масштабуванням швидкості

Залишена сама на себе, система МД зберігає повну енергію (мікроканонічний ансамбль, NVE), але температура дрейфує, поки кінетична й потенціальна енергія обмінюються. Щоб утримати цільову температуру (канонічний ансамбль, NVT), найпростіший термостат масштабує кожну швидкість спільним множником щокроку, щоб миттєва кінетична температура точно збігалась з цільовою:

function rescaleVelocities(particles, targetT) {
  const n = particles.length;
  const kineticSum = particles.reduce((s, p) => s + p.vx**2 + p.vy**2, 0);
  const currentT = kineticSum / n;
  if (currentT < 1e-9) return;
  const lambda = Math.sqrt(targetT / currentT);
  for (const p of particles) { p.vx *= lambda; p.vy *= lambda; }
}

8. Повний цикл симуляції

Складання обчислення сил, інтегрування, меж і термостатування в одну функцію на кадр:

function computeForces(particles, boxSize, cutoff = 2.5, epsilon = 1, sigma = 1) {
  for (const p of particles) { p.ax = 0; p.ay = 0; }
  for (let i = 0; i < particles.length; i++) {
    for (let j = i + 1; j < particles.length; j++) {
      let dx = particles[j].x - particles[i].x;
      let dy = particles[j].y - particles[i].y;
      ({ dx, dy } = minimumImageDelta(dx, dy, boxSize));
      const r2 = dx * dx + dy * dy;
      if (r2 > cutoff * cutoff || r2 < 1e-6) continue;
      const r = Math.sqrt(r2);
      const fMag = ljForceMagnitude(r, epsilon, sigma) / r; // заздалегідь ділимо, щоб спроєктувати на dx, dy
      particles[i].ax -= fMag * dx; particles[i].ay -= fMag * dy;
      particles[j].ax += fMag * dx; particles[j].ay += fMag * dy;
    }
  }
}

function simulate(particles, boxSize, dt, targetT) {
  velocityVerletStep(particles, null, dt, p => computeForces(p, boxSize));
  for (const p of particles) wrapPosition(p, boxSize);
  if (targetT !== null) rescaleVelocities(particles, targetT);
}

Запустіть цей цикл при низькій густині й високій T* — отримаєте невпорядкований газ; знизьте температуру або підвищте густину — і частинки конденсуються в краплі рідини, а потім у кристалічне тверде тіло — та сама фазова поведінка, розглянута якісно у статті сайту про рівняння Ван-дер-Ваальса і кількісно на фазовій діаграмі цієї системи ЛД.

9. Типові пастки

• Розбіжність сили при r→0: член r⁻¹³ у силі ЛД злітає в нескінченність, якщо дві частинки колись точно перекриються — завжди стартуйте з ґратки без перекриття і використовуйте достатньо малий крок часу, щоб частинки ніколи не зближувались менш ніж на ≈0.8σ.
• Забути найближче зображення: обчислення прямого dx = x_j − x_i без загортання змушує частинки біля протилежних країв ящика відчувати величезну, неправильну далекодіючу силу замість коректної короткодіючої.
• Занадто великий крок часу: крутий відштовхувальний кор ЛД потребує малого dt (типово 0.001–0.005 у зведених одиницях часу) — дрейф енергії вгору протягом прогону — класичний симптом завеликого кроку.
• Пропуск зсуву потенціалу: необрізаний зсув вводить малу розривність сили точно при r = r_c, що повільно вносить енергію протягом тисяч кроків.
• Масштабування щоразового кроку: постійне жорстке масштабування пригнічує природні флуктуації температури; масштабуйте лише кожні N кроків або на початковій фазі зрівноважування, а потім вимикайте (або переходьте на м'якший термостат), щоб семплювати коректний ансамбль.