Динаміка Ланжевена та термостати

1. Навіщо молекулярній динаміці термостат

Рівняння руху Ньютона, інтегровані симплектичною схемою на кшталт швидкісного Верле, зберігають повну енергію з високою точністю. Це дає мікроканонічний ансамбль (NVE): фіксовані кількість частинок N, об'єм V та енергія E. Це природний результат чистої ньютонівської динаміки, але рідко саме те, що потрібно для порівняння з реальним експериментом, який майже завжди проводиться за сталої температури, а не сталої енергії — пробірка у водяній бані, білок у розчині, газ у кімнаті.

Канонічний ансамбль (NVT) натомість фіксує N, V та T, дозволяючи енергії коливатись, поки система обмінюється теплом з неявним резервуаром. Термостат — це частина симуляційного рушія, яка відтворює цей обмін без буквальної симуляції мільярдів молекул розчинника чи зіткнень зі стінками. Два домінантні сімейства:

2. Рівняння Ланжевена: тертя і шум

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

Фізично γ представляє швидкість, з якою неявна баня забирає кінетичну енергію через зіткнення, тоді як ξ(t) представляє випадкові поштовхи імпульсу від тих самих зіткнень. Самі по собі, будь-який з доданків штовхав би систему від рівноваги — чисте тертя заморозило б увесь рух до v = 0, чистий шум нагрівав би систему безмежно. Лише правильна комбінація обох дає стабільну рівноважну температуру.

3. Флуктуаційно-дисипаційне обмеження

Коефіцієнт тертя γ і амплітуду шуму не можна обирати незалежно, якщо потрібно, щоб система врівноважувалась за конкретної цільової температури T. Вони жорстко пов'язані тим самим співвідношенням, що й у флуктуаційно-дисипаційній теоремі:

Саме тому динаміка Ланжевена є коректним термостатом, а не довільним хаком: саме тертя порушило б детальну рівновагу і забрало б всю енергію із системи; додавання шуму точно з силою флуктуаційно-дисипаційного співвідношення відновлює детальну рівновагу відносно канонічного (больцманівського) розподілу, тож усереднення в часі вздовж довгої траєкторії Ланжевена дорівнюють середнім по канонічному ансамблю.

4. Вибір коефіцієнта тертя

γ — вільний параметр, який контролює компроміс між силою термостатування і відповідністю базовій ньютонівській динаміці:

• Мале γ (слабкий зв'язок): траєкторія залишається близькою до справжньої NVE-траєкторії між рідкісними корекціями термостата, добре зберігаючи динамічні величини на кшталт коефіцієнтів дифузії та автокореляційних функцій швидкості, але дрейф температури коригується повільно.
• Велике γ (сильний зв'язок): температура контролюється жорстко, а рівновага досягається швидко, але динаміка стає надмірно демпфованою — ближчою до броунівського руху, ніж до балістичного молекулярного руху — що спотворює будь-яку динамічну властивість, яку хочеться виміряти.
• Практичне правило: обирати 1/γ у кілька разів більшим за найкоротший час релаксації властивості, яку важливо зберегти; для суто статичних/структурних величин (радіальні функції розподілу, вільні енергії) велике γ цілком підходить і часто навіть краще для швидшого врівноваження.

5. Термостат Нозе-Гувера

Нозе (1984) і Гувер (1985) обрали інший підхід: замість вприскування випадкового шуму, додати одну додаткову детерміновану змінну, що представляє ефективну силу зв'язку з тепловою банею, і дозволити її власному рівнянню руху штовхати кінетичну енергію системи до цільового значення:

Коли миттєва кінетична енергія K перевищує ціль 3Nk_BT/2, ξ зростає і додаткове тертя сповільнює частинки; коли K падає нижче цілі, ξ може стати від'ємним і фактично прискорити частинки. Цей контур зворотного зв'язку повністю детермінований і оборотний у часі, і — головне — можна показати, що він семплує точний канонічний розподіл, а не просто контролює середню температуру.

Параметр Q задає власну інерцію термостата: мале Q дає "жорсткий", швидкореагуючий термостат, здатний перерегулювати і осцилювати; велике Q дає "м'який" термостат, що повільно дрейфує до цільової температури, ближчий у короткостроковій перспективі до мікроканонічної динаміки.

6. Ланцюги Нозе-Гувера та ергодичність

Одна змінна Нозе-Гувера має відомий режим відмови: для малих чи жорстких систем (класичний приклад — одиночний гармонічний осцилятор) поєднана динаміка системи-з-термостатом може стати квазіперіодичною, а не ергодичною, осцилюючи навколо цільової температури, так і не досягнувши правильного семплування повного канонічного розподілу.

Виправлення Мартини, Кляйна і Такермана — ланцюг Нозе-Гувера — з'єднує початкову термостатну змінну ξ₁ з другою термостатною змінною ξ₂, яка гасить власні флуктуації ξ₁, і яку своєю чергою можна з'єднати з ξ₃, і так для короткого ланцюга (типова довжина 3-5). Кожна ланка поглинає залишкову неергодичність попередньої, і на практиці короткий ланцюг відновлює правильне канонічне семплування практично для будь-якого розміру системи, за незначної додаткової обчислювальної вартості.

7. JavaScript-інтегратор (BAOAB)

// Схема розщеплення BAOAB для динаміки Ланжевена — семплує точний
// канонічний розподіл навіть при досить великих кроках за часом.
function baoabStep(x, v, m, force, gamma, kT, dt, randn) {
  const halfDt = dt / 2;

  // B: напівкрок імпульсу від сил
  v += (force(x) / m) * halfDt;

  // A: напівкрок дрейфу позиції
  x += v * halfDt;

  // O: точне оновлення Орнштейна-Уленбека — тертя і шум разом
  const c1 = Math.exp(-gamma * dt);
  const c2 = Math.sqrt((1 - c1 * c1) * kT / m);
  v = c1 * v + c2 * randn();

  // A: другий напівкрок дрейфу
  x += v * halfDt;

  // B: другий напівкрок від (оновлених) сил
  v += (force(x) / m) * halfDt;

  return { x, v };
}

// Нозе-Гувер: детермінований зворотний зв'язок, без випадковості
function noseHooverStep(x, v, xi, m, force, Q, kT, dof, dt) {
  const K = 0.5 * m * v * v;
  const dxi = (2 * K - dof * kT) / Q;
  xi += dxi * dt;
  v += (force(x) / m - xi * v) * dt;
  x += v * dt;
  return { x, v, xi };
}

8. Застосування та підводні камені

• Симуляція білків і біомолекул: термостати Ланжевена з помірним тертям — типовий вибір для врівноваження сольватованих біомолекул, бо вони надійні й дешеві; продуктивні прогони для обчислення динамічних властивостей часто переходять на ланцюги Нозе-Гувера або слабший зв'язок.
• Транспортні коефіцієнти: обидва термостати можуть спотворювати транспортні коефіцієнти Ґріна-Кубо за надто сильного зв'язку — завжди перевіряйте, чи результати збігаються при послабленні зв'язку термостата.
• «Літаючий крижаний кубик»: відомий артефакт Нозе-Гувера, коли енергія нефізично перетікає в трансляційні/обертальні моди центру мас майже жорсткої молекули, вимагаючи періодичного видалення руху центру мас як окремої корекції.
• Малі системи: один термостат Нозе-Гувера може не бути ергодичним для жорстких, малорозмірних систем — використовуйте ланцюг Нозе-Гувера (див. розділ 6) або перейдіть на Ланжевен.
• Зв'язок з Монте-Карло: симуляція моделі Ізінга на сайті семплує той самий канонічний (больцманівський) розподіл дискретними кроками Метрополіса замість неперервної динаміки — корисний контраст, як досягти однієї й тієї самої статистичної цілі зовсім іншим шляхом.