Молекулярна динаміка: потенціал Леннарда-Джонса

3D Simulations

Хімія · Статистична механіка · Симуляція

📅 Липень 2026 ⏱ ≈ 10 хв читання 🎯 Середній рівень · Оновлено: 9 липня 2026

Молекулярна динаміка: потенціал Леннарда-Джонса

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

Автор: MySimulator Team · Рецензія: Редакційна рецензія MySimulator

1. Потенціал Леннарда-Джонса

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

U(r) = 4ε \cdot [ (σ/r)¹² - (σ/r)⁶ ] ε = глибина потенціальної ями (сила зв'язку, в одиницях енергії) σ = відстань, на якій U(r) = 0 (приблизно діаметр атома) r = відстань між центрами двох частинок

Притягальний доданок −(σ/r)⁶ обраний не довільно — це справжній результат провідного порядку теорії збурень квантової механіки для взаємодії наведений-диполь/наведений-диполь між двома поляризовними атомами. Відштовхувальний доданок (σ/r)¹² не має настільки строгого виведення; дванадцяту степінь обрали здебільшого через дешевизну обчислення (це просто квадрат шостого степеня) та тому, що вона задовільно відтворює експериментальні дані стисливості.

Потенціал має мінімум при r_min = 2^1/6·σ, де U(r_min) = −ε. Це природна рівноважна відстань між двома атомами за відсутності інших сил — точка, де притягання та відштовхування точно врівноважуються.

Типові параметри: для аргону ε/k_B ≈ 120 К (ε ≈ 1,65×10⁻²¹ Дж) та σ ≈ 3,4 Å. Цих двох чисел достатньо, щоб відтворити температуру плавлення, кипіння та критичну точку аргону з точністю до кількох відсотків — вражаюча віддача від моделі з двома параметрами.

2. Як з потенціалу отримати силу

Молекулярна динаміка інтегрує другий закон Ньютона, тож нам потрібна не потенціальна енергія, а сила. Сила між парою частинок — це від'ємний градієнт U за відстанню:

F(r) = -dU/dr = 24ε/r \cdot [ 2(σ/r)¹² - (σ/r)⁶ ] F(r) > 0 \to відштовхування (розштовхує атоми на коротких відстанях) F(r) < 0 \to притягання (стягує атоми на великих відстанях)

У трьох вимірах вектор сили на частинку i від частинки j спрямований уздовж вектора розділення r_ij = r_i − r_j:

F⃗ᵢⱼ = F(r) \cdot (r⃗ᵢⱼ / r) (одиничний вектор уздовж зв'язку) F⃗ⱼᵢ = -F⃗ᵢⱼ (третій закон Ньютона)

Для N частинок повна сила на частинку i — це попарна сума за всіма іншими частинками:

F⃗ᵢ = Σⱼ\neqᵢ F⃗ᵢⱼ

Такий наївний підсумок коштує O(N²) попарних обчислень на крок — прийнятно для кількох сотень частинок у браузерній демонстрації, але саме тому промислові пакети (LAMMPS, GROMACS) використовують списки клітинок та списки сусідів, щоб знизити вартість приблизно до O(N) після введення радіуса обрізання (див. розділ 4).

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

Знаючи сили, потрібно просувати положення та швидкості в часі. Робочий алгоритм молекулярної динаміки — швидкісний Верле (velocity Verlet), цінний тим, що він оборотний у часі, симплектичний (зберігає обсяг фазового простору й, з дуже гарним наближенням, повну енергію на довгих прогонах) і дешевий — одне обчислення сили на крок:

r⃗(t+Δt) = r⃗(t) + v⃗(t)\cdotΔt + ½\cdota⃗(t)\cdotΔt² a⃗(t+Δt) = F⃗(r⃗(t+Δt)) / m (перерахувати сили в нових положеннях) v⃗(t+Δt) = v⃗(t) + ½\cdot[a⃗(t) + a⃗(t+Δt)]\cdotΔt

Порівняйте це зі звичайним інтегруванням Ейлера (v += a·Δt; r += v·Δt), яке не є оборотним у часі й систематично зносить енергію вгору або вниз протягом тисяч кроків — воно помітно "нагріє" або "заморозить" LJ-рідину, яка мала б лишатися при сталій енергії. Локальна похибка швидкісного Верле — O(Δt⁴) на крок (глобальна O(Δt²)), порівняно з O(Δt²) локально для Ейлера, за практично однакової обчислювальної вартості.

Вибір кроку часу: Δt має розв'язувати найшвидше коливання в системі — для LJ-рідини поблизу мінімуму потенціалу це приблизно період коливання на дні ями. Загальне правило: Δt ≈ 0,005–0,01 у зведених LJ-одиницях (τ = σ√(m/ε)), тобто близько 1–2 фемтосекунд для повноатомних симуляцій з явними атомами водню.

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

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

якщо x > L: x -= L якщо x < 0: x += L (циклічне обгортання в кожному вимірі)

Оскільки U(r) спадає як r⁻⁶, її внесок за межами кількох σ нехтовний. Промислові пакети застосовують радіус обрізання r_c (зазвичай 2,5σ), за яким взаємодію просто прирівнюють до нуля, разом із конвенцією найближчого зображення: для кожної пари враховується лише найближча періодична копія частинки j при обчисленні r_ij.

Обрізання + списки сусідів: список Верле кешує, які частинки перебувають у межах r_c плюс "запасна" відстань, перебудовується лише раз на кілька кроків — знижуючи цикл обчислення сил приблизно до O(N) за фіксованої густини.
Списки клітинок: бокс поділяється на клітинки розміром ≥ r_c; частинці треба перевіряти лише власну клітинку та 26 сусідніх (у 3D), даючи ту саму асимптотику з простішою структурою даних.
Поправки на "хвіст": оскільки обрізання відкидає невеликий притягальний внесок, до тиску та енергії зазвичай додають аналітичну поправку на далекий "хвіст", припускаючи однорідну густину за межами r_c.

Симуляція молекулярної динаміки

Спостерігайте, як LJ-газ конденсується в краплю рідини або кристалізується в решітку, змінюючи густину та температуру в реальному часі.

Відкрити симуляцію →

5. Термостати та контроль температури

Швидкісний Верле сам по собі зберігає повну енергію (мікроканонічний ансамбль, NVE). Більшість фізичних питань — що відбувається за 300 К? у точці плавлення? — вимагають натомість фіксувати температуру (канонічний ансамбль, NVT), а для цього потрібен термостат.

Миттєва температура визначається через теорему про рівнорозподіл із кінетичної енергії всіх N частинок у 3D:

T = (1/3Nk_B) \cdot Σᵢ mᵢ\cdotvᵢ²

Масштабування швидкостей: найгрубіший термостат — щокроку множити кожну швидкість на √(T_ціль/T_{поточна}). Простий, але не відтворює канонічний ансамбль коректно (немає реалістичних флуктуацій).
Термостат Берендсена: плавно масштабує швидкості до цілі з часом релаксації τ, уникаючи різких стрибків миттєвого масштабування, але все ще не строго канонічний (надто сильно пригнічує флуктуації температури — проблема "летючого крижаного кубика" у крайніх випадках).
Термостат Нозе-Гувера: вводить додаткову динамічну змінну, що представляє ступінь свободи теплової бані, зв'язану з фізичними частинками через розширений лагранжіан. Це дає математично коректний канонічний (NVT) розподіл і є стандартним вибором у промислових пакетах.
Термостат Ланжевена: додає тертя та випадкову (стохастичну) силу до кожної частинки, імітуючи зіткнення з неявним розчинником чи тепловою банею. Простий у реалізації та фізично інтуїтивний, ціною зміни динаміки (він не зберігає імпульс так, як Нозе-Гувер).

6. Чому це досі важливо

Благородні гази та прості рідини

Потенціал LJ залишається підручниковою моделлю для благородних газів (Ar, Kr, Xe) та стандартним еталоном для перевірки нових інтеграторів, термостатів і паралельних алгоритмів МД — саме тому, що його фазова діаграма (потрійна точка, критична точка, крива плавлення) відома з високою точністю завдяки десятиліттям симуляційних досліджень.

Силові поля для біомолекул

Сучасні повноатомні силові поля, як-от AMBER, CHARMM та OPLS, досі використовують член LJ для кожної незв'язаної пари атомів у поєднанні з кулонівською електростатикою та доданками зв'язків/кутів/торсій. Форма 12-6 зберігається здебільшого з історичних та обчислювальних причин — заміна на фізично точніше експоненційне відштовхування (потенціал Бекінгема) коштовніша при кожному обчисленні й рідко змінює результати настільки, щоб виправдати ціну в масштабі.

Матеріалознавство

LJ та його варіанти моделюють тверді благороднагазові тіла, деякі прості метали (через розширення методом вбудованого атома), і часто використовуються як спрощений proxy-потенціал у розробці методів для дослідження нуклеації, руйнування та меж зерен перед переходом до дорожчих квантових чи багаточастинкових потенціалів.

Часті запитання

Чому коди МД використовують степеневий закон 12-6, а не точніше відштовхування?

Відштовхування r⁻¹² зручне для обчислень (це просто квадрат доданку r⁻⁶, тож обидва можуть спільно використовувати те саме обчислення відстані), і воно достатньо добре відтворює об'ємні термодинамічні властивості простих рідин для більшості цілей. Точніші експоненційні (типу Бекінгема) відштовхування існують і застосовуються, коли важлива точність на дуже коротких відстанях, але вони коштують додаткових трансцендентних обчислень на пару, що суттєво накопичується при мільярдах парних взаємодій у великій симуляції.

Що трапиться, якщо крок часу занадто великий?

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

Як обчислюється тиск у симуляції Леннарда-Джонса?

Тиск обчислюється через теорему віріала, яка пов'язує ансамблево-усереднену кінетичну енергію та суму скалярних добутків попарних сил на розділення за всіма парами частинок: P·V = N·k_B·T + (1/3)·⟨Σᵢ<ⱼ F⃗ᵢⱼ·r⃗ᵢⱼ⟩. Перший доданок — це внесок ідеального газу; другий, віріальний доданок, охоплює поправку від міжчастинкових сил і може бути додатним (чисте відштовхування, за високої густини) або від'ємним (чисте притягання, за низької густини поблизу критичної точки).

Хто такий Джон Леннард-Джонс і коли він запропонував цей потенціал?

Сер Джон Едвард Леннард-Джонс — британський математик і теоретичний хімік, який запропонував потенціал 12-6 у 1924 році, вивчаючи рівняння стану газів у Брістолі, а пізніше став першим професором теоретичної хімії в Кембриджі. Його оригінальна робота використовувала потенціал для узгодження другого віріального коефіцієнта газів, як-от аргон, за десятиліття до перших цифрових симуляцій МД (що почалися у 1950-х із моделей твердих сфер і прийняли взаємодії LJ у 1960-х, зокрема в знаковій симуляції рідкого аргону Рахмана 1964 року).

Чим молекулярна динаміка відрізняється від симуляції Монте-Карло?

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

Чому LJ-система кристалізується за низької температури?

За низької температури кінетична енергія на частинку (порядку k_BT) мала порівняно з глибиною потенціальної ями ε. Частинки осідають поблизу мінімуму енергії r_min = 2^(1/6)σ своїх сусідів, а найщільніший спосіб упакувати сфери за фіксованої низької енергії — це щільно упакована решітка (ГЦК або ГЩУ). Із підвищенням температури тепловий рух долає глибину ями, решітка плавиться в рідину, а зрештою (за низької густини) система стає розрідженим газом — та сама структура потрійної точки, що спостерігається в реальних благородних газах.

Чи може потенціал Леннарда-Джонса моделювати хімічні зв'язки?

Не безпосередньо. LJ описує незв'язані, ван-дер-ваальсові взаємодії між атомами, які не пов'язані один з одним хімічно. Ковалентні зв'язки всередині молекули натомість моделюються гармонічними (або морзевими) потенціалами розтягу зв'язку, доданками згинання кутів та торсійними доданками; член LJ залишають для взаємодій між атомами в різних молекулах або між атомами тієї самої молекули, розділеними трьома чи більше зв'язками (так звані взаємодії 1-4 та далі).

Які одиниці зазвичай використовують у LJ-симуляціях?

Щоб уникнути роботи з крихітними СІ-числами, коди МД зазвичай використовують "зведені" LJ-одиниці, де довжина вимірюється в одиницях σ, енергія — в одиницях ε, а маса — в одиницях маси частинки m. Час тоді природно вимірюється в одиницях τ = σ√(m/ε), а температура — в одиницях ε/k_B. У цих одиницях ε = σ = m = 1, і всі величини стають безрозмірними числами порядку одиниці, що зручно для обчислень і робить результати безпосередньо переносними між різними речовинами простим масштабуванням.