Практикум вычислительные методы исследования молекулярной динамики аксенова е. В., Кшевецкий м. С. Учебнометодическое пособие (для студентов физического факультета)
 Скачать 466.62 Kb.
|
1 2 выбираются ли частицы случайно или последовательно? 37 7. Моделирование микроканонического ансамбля методом Монте-Карло проводится при фиксированной полной энергии независимо от темпе- ратуры. Определим температуру соотношением (14). Используя это соотношение, получите температуру. Связана ли полученная темпера- тура со средней энергией демона? 8. Добавьте к программе подпрограмму, которая вычисляет распределе- ние вероятностей P (E d ) демона и его среднюю энергию. Нарисуйте график логарифма P (E d ) как функцию от E d и проверьте формулу (27). Оцените соответствующую величину k B T . 9. Определите величину k B T из выражения (43). Совпадают ли ваши две оценки k B T ? Лабораторная работа 2.3. Двумерная модель Изинга Целью настоящей работы является исследование двумерной модели Изинга методом Монте-Карло в каноническом ансамбле. Содержание работы 1. Выпишите математическую модель (выражения для энергии в модели Изинга) для моделирования методом Монте-Карло. Определите набор входных параметров. 2. Выполните обезразмеривание уравнений и начальных параметров. 3. Выполните программу. Найдите средние значения энергии и намагни- ченности. Здесь и в других заданиях расчеты проводите путем усред- нения по последним k шагам метода Монте-Карло (самостоятельно модифицируйте программу). Значения k выбирайте самостоятельно таким образом, чтобы средние значения от серии к серии изменялись несильно. 4. Найдите среднее значение поля, при котором в среднем половина ча- стиц имеет спин, смотрящий против поля. Найдите при заданном поле среднее число частиц со спинами, направленными вверх и вниз. 38 5. Задания 3 и 4 выполните при разных J (положительных и отрица- тельных) при одной и той же температуре. 6. Реализуйте в программе возможность изменения температуры и смо- делируйте фазовый переход. Постройте зависимости средних энергии, намагниченности, теплоемкости и восприимчивости от температуры. Приложение. Комментарии к коду программ Прилагаемые программы monte-carlo-metropolis.c для работы 2.1. и monte-carlo-isolate.c для работы 2.2. на языке С снабжены подробны- ми комментариями. Программы реализуют моделирование методом Монте- Карло, соответственно, канонического и микроканонического ансамблей в безразмерных переменных. Сначала приведены функции, вычисляющие энергию взаимодействия одной частицы, кинетическую и потенциальную энергию системы. Пред- полагается, что на систему наложены периодические граничные условия. При расчете энергии взаимодействия учитывается взаимодействие только с ближайшей копией каждой из частиц. Особое внимание следует обратить на функции metropolis_step в ра- боте 2.1. и monte_carlo_step в работе 2.2., реализующие расчет координат и скоростей для одного шага метода Монте-Карло. Здесь параметр delta задает, на какую долю можно изменять координаты и скорости частицы. В этой функции принимается решение об изменении или неизменении со- стояния системы в зависимости от энергии демона и энергии пробного из- менения состояния системы. Далее следуют функции, которые собирают данные для построения рас- пределения по скоростям и осуществляют запись и чтение данных о коор- динатах и скоростях в файл. Эти функции можно не изучать подробно. Важная функция initialize задает начальные координаты и скорости частиц. В этой функции следует разобраться подробно. По координатам ча- стицы равномерно распределяются по объему системы (разберитесь, как), а проекции скоростей задаются случайным образом в пределах от −Vmax до +Vmax. Также учитывается условие, что суммарный импульс системы 39 равен нулю. Далее идет вспомогательный блок, необходимый для графического отоб- ражения распределения по координатам и скоростям. Наконец, основной является функция main, где собственно и происходит выполнение моделирования. При запуске программы следует ввести число частиц N , максимальное значение проекции скорости v max и допустимую долю для расчета измене- ний координат и скорости частиц δ. В случае работы 2.1. также потребует- ся ввести температуру T . Начнется процесс моделирования. При нажатии комбинации клавиш Ctrl+C, вычисления приостановятся. Далее следуйте инструкциям на экране. Прилагаемая программа ising-metropolis.c для работы 2.3. на языке С снабжена подробными комментариями. Программа реализуют модели- рование двумерной модели Изинга методом Монте-Карло в безразмерных переменных. Сначала приведены функции, вычисляющие сумму спинов ближайших соседей, полную энергию системы и намагниченность. Предполагается, что на систему наложены периодические граничные условия. Особое внимание следует обратить на функцию metropolis_step, реа- лизующую расчет спинов частиц и энергии системы для одного шага ме- тода Монте-Карло. В этой функции принимается решение об изменении или неизменении состояния системы в зависимости от энергии пробного изменения состояния системы. Далее следуют функции, которые собирают данные для построения спи- нов частиц и осуществляют запись и чтение данных о спинах в файл. Эти функции можно не изучать подробно. Важная функция initialize задает начальные значения спинов ча- стиц. Спины задаются случайным образом. Далее идет вспомогательный блок, необходимый для графического отоб- ражения распределения спинов частиц. Наконец, основной является функция main, где собственно и происходит выполнение моделирования. При запуске программы следует ввести размеры ячейки по осям x и 40 y, константу обменного взаимодействия J , внешнее поле H и температу- ру T . Начнется процесс моделирования. При нажатии комбинации кла- виш Ctrl+C, вычисления приостановятся. Далее следуйте инструкциям на экране. 3 Броуновская динамика 3.1 Уравнение Ланжевена При рассмотрении броуновского движения взвешенной частицы исполь- зуется метод Ланжевена, в котором уравнение движения записывается в виде m d 2 r dt 2 = − 1 B dr dt + F(t) (44) или m dv dt = − v B + F(t), (45) где m — масса броуновской частицы, v/B — сила трения, пропорциональ- ная скорости, F(t) — случайная (стохастическая) сила, обладающая рядом статистических свойств. В простейшем случае (приближение белого шума) на случайную силу F(t) накладывают следующие ограничения hF α (t)i = 0, (46) hF α (t 1 )F β (t 2 )i = Cδ(t 1 − t 2 )δ αβ (47) Здесь угловыми скобками обозначено усреднение по ансамблю реализаций функции F(t), индексы α и β нумеруют проекции случайной силы, δ(t 1 − t 2 ) — дельта-функция Дирака, а δ αβ — символ Кронекера. Функция hF α (t 1 )F β (t 2 )i определяет степень статистической независимо- сти величин F α (t 1 ) и F β (t 2 ) и называется корреляционной функцией. С учетом эргодической гипотезы усреднение по ансамблю может быть за- менено на усреднение по времени, тогда корреляционная функция может быть представлена в виде hF α (t + τ )F β (t)i = lim t→∞ 1 t Z t 0 F α (t 0 + τ )F β (t 0 )dt 0 (48) 41 Возвращаясь назад к соотношениям (46) и (47), видим, что они описыва- ют ситуацию, когда действия случайной силы в разные моменты времени статистически независимы, а сама случайная сила в среднем равна нулю. 3.2 Решение уравнения Ланжевена Рассмотрим далее проекцию уравнения Ланжевена на ось x системы координат. Построим сначала решение однородного уравнения (45) c F = 0. Получим v x (t) = u exp − t mB (49) Неоднородное уравнение c F 6= 0 решим, считая, что u = u(t). В качестве начального условия при t = 0 возьмем v x = v 0x . Имеем v x (t) = v 0x exp − t mB + 1 m Z t 0 F x (τ ) exp τ − t mB dτ. (50) Усредним по ансамблю реализаций случайных сил F x (t). Учитывая свой- ство (46), находим hv x (t)i = v 0x exp − t mB + 1 m Z t 0 hF x (τ )i exp τ − t mB dτ = = v 0x exp − t mB . (51) Теперь рассмотрим hv 2 x (t)i. Возведем (50) в квадрат и запишем произ- ведение интегралов как двойной интеграл. Получаем hv 2 x (t)i = v 2 0x exp − 2t mB + 2 v 0x m exp − 2t mB Z t 0 hF x (τ )i exp τ mB dτ + + 1 m 2 exp − 2t mB Z t 0 Z t 0 hF x (τ )F x (τ 1 )i exp τ + τ 1 mB dτ 1 dτ. (52) Первый интеграл в (52) равен нулю (см. (46)), а двойной интеграл легко вычисляется при использовании свойства (47) корреляционной функции: C Z t 0 Z t 0 δ(τ − τ 1 ) exp τ + τ 1 mB dτ 1 dτ = C mB 2 exp 2t mB − 1 (53) 42 В результате hv 2 x (t)i = v 2 0x exp − 2t mB + CB 2m 1 − exp − 2t mB (54) При достаточно большом времени (t → ∞) влияние начальных условий сглаживается. При этом hv 2 x (t)i принимает равновесное значение, равное k B T /m. Поэтому с помощью (54) находим hv 2 x (t)i → hv 2 x i = CB 2m = k B T m (55) Отсюда получаем C = 2k B T B Из выражения (54) следует, что равновесное значение среднего квадрата скорости броуновской частицы установится лишь по прошествии времени, значительно превышающего t mB. Найдем теперь зависимость координаты броуновской частицы от време- ни. Воспользуемся соотношением (50) для скорости броуновской частицы, получаемом при интегрировании уравнения (45). Интегрируя выражение (50) для v x (t) по времени при начальном условии x(0) = x 0 , получаем x(t) = x 0 + Z t 0 v x (τ )dτ (56) или x(t) = x 0 + mBv 0x 1 − exp − t mB + + 1 m Z t 0 Z τ 1 0 F x (τ ) exp τ − τ 1 mB dτ dτ 1 . (57) Меняя местами порядок интегрирования и снимая один интеграл, получаем x(t) = x 0 + mBv 0x 1 − exp − t mB + B Z t 0 F x (τ ) 1 − exp τ − t mB dτ. (58) Возведем равенство (58) в квадрат и усредним по ансамблю реализаций случайных сил F x (t). Получим hx 2 (t)i = x 0 + mBv 0x 1 − exp − t mB 2 + + B 2 C Z t 0 Z t 0 δ(τ − τ 1 ) 1 − exp τ − t mB 1 − exp τ 1 − t mB dτ dτ 1 . (59) 43 Снимем один интеграл за счет дельта-функции и вычислим оставшийся интеграл. Получаем hx 2 (t)i = = x 2 0 + 2x 0 v 0x mB 1 − exp − t mB + v 2 0x m 2 B 2 1 − exp − t mB 2 + + CB 3 m 2 4 exp − t mB − exp − 2t mB − 3 + CB 2 t. (60) 3.3 Численное решение уравнения Ланжевена Численное решение уравнения Ланжевена можно построить следующим образом. Введем шаг по времени ∆t. В начальный момент времени поло- жим t = 0. Также будем считать, что координаты и проекции скорости частицы в начальный момент времени равны нулю, то есть, r 0 = 0, v 0 = 0. (61) Нас будут интересовать отсчеты координат r i и скоростей v i в моменты времени t i = i∆t, где i = 1, . . . , n, а n — количество шагов по времени. Для решения уравнения (44) мы будем использовать метод Эйлера пер- вого порядка точности. В этом случае численная схема для нахождения координат и скоростей броуновской частицы в заданные моменты времени имеет вид r i+1 = r i + v i ∆t, v i+1 = v i 1 − ∆t mB + F i ∆t m , (62) где вектор F i представляет собой случайную силу действующую в момент времени t i . Мы будем считать, что каждая из компонент этой силы пред- ставляет собой случайную величину равномерно распределенную в диапа- зоне (−F max , F max ). Рекуррентная схема (62) c начальными условиями (61) позволяет нам выполнить один проход моделирования и найти координаты и скорости броуновской частицы для одной единственной реализации случайной силы F(t). Однако для изучения статистических свойств броуновской частицы 44 этого недостаточно, поэтому мы будем повторять вычисления по схеме (62) c начальными условиями (61) многократно и накапливать полученные ре- зультаты. Так, например, для того, чтобы найти средние значения квадрата координаты и скорости броуновской частицы в каждый момент времени, мы будем суммировать квадраты координат и скоростей частицы в каждый момент времени для всех проходов S M r 2 i = M X j=1 (r (j) i ) 2 , S M v 2 i = M X j=1 (v (j) i ) 2 , (63) где r (j) i и v (j) i — значения координат и скорости частицы в момент времени t i на проходе j, а M — количество выполненных проходов. После выпол- нения каждого прохода мы будем оценивать относительную погрешность вычисления средних величин с помощью следующих формул ∆ r = max i S M r 2 i − M · (r (M ) i ) 2 (M − 1)S M r 2 i , ∆ v = max i S M v 2 i − M · (v (M ) i ) 2 (M − 1)S M v 2 i (64) Если наибольшее из значений ∆ r и ∆ v окажется меньше наперед заданного значения ε, то мы достигли требуемой точности расчетов, и можно далее не увеличивать количество проходов M . 3.4 Анализ особенностей реализации случайной силы в численной схеме В предыдущем параграфе мы ввели силу F(t), действие которой в тече- нии интервала времени [t i , t i+1 ) определяется случайной величиной F i . В этом параграфе мы будем изучать статистические свойства этой случайной силы. Изучим сперва статистические свойства проекций случайной величины F i , а затем вернемся к изучению статистических свойств, задаваемой ими случайной силы F(t). Мы начнем с вычисления среднего значения проек- ции F x,i случайной величины F i на ось x (проекции на другие оси координат находятся аналогично). Вспоминая, что F x,i представляет собой случайную величину ξ равномерно распределенную в интервале (−F max , F max ), имеем ρ(ξ) = 1 2F max , (65) 45 где ρ(ξ) — плотность вероятности распределения случайной величины ξ. Тогда hF α,i i = hF x,i i = Z +F max −F max ξρ(ξ)dξ = 0. (66) Теперь изучим парный коррелятор hF α,i F β,j i. Для начала заметим что случайная величина F i не зависит от случайной величины F j при i 6= j, кроме того мы считаем, что проекции на разные оси случайной величины F i статистически независимы, тогда hF α,i F β,j i = 0, если i 6= j или α 6= β. (67) Итого нам осталось изучить всего один случай, когда i = j и α = β. Мы будем изучать только hF 2 x,i i, а для других направлений результат будет тем же самым, итак hF 2 α,i i = hF 2 x,i i = Z +F max −F max ξ 2 ρ(ξ)dξ = F 2 max 3 (68) Объединяя (67) и (68) получаем финальное выражение для парного кор- релятора hF α,i F β,j i = F 2 max 3 δ αβ δ i,j (69) Теперь мы можем вернуться к изучению случайной силы F(t). По опре- делению имеем F(t) = ∞ X i=0 F i θ(t − t i ) θ(t i+1 − t), (70) где θ(·) представляет собой θ-функцию Хевисайда. Тогда hF(t)i = * ∞ X i=0 F i θ(t − t i ) θ(t i+1 − t) + = = ∞ X i=0 hF i i θ(t − t i ) θ(t i+1 − t) = 0. (71) 46 Рассмотрим теперь парную корреляционную функцию hF(t)F(t + τ )i = = * ∞ X i=0 ∞ X j=0 F i F j θ(t − t i )θ(t i+1 − t)θ(t + τ − t i )θ(t i+1 − t − τ ) + = = d F 2 max 3 * ∞ X i=0 θ(t − t i )θ(t i+1 − t)θ(t + τ − t i )θ(t i+1 − t − τ ) + , (72) где d — размерность пространства. Несложно увидеть, что выражение в угловых скобках отлично от нуля только если |τ | < ∆t. Переходя от усред- нения по ансамблю к усреднению по времени и учитывая периодичность выражения в угловых скобках, получаем hF(t)F(t + τ )i = d F 2 max 3 × 1 + τ /∆t, если − ∆t < τ ≤ 0; 1 − τ /∆t, если 0 < τ < ∆t; 0, если |τ | ≥ ∆t. (73) 3.5 Сила, действующая продолжительное время До сих пор мы рассматривали ситуацию, когда сила, действующая на ча- стицу, была случайной и действовала только в один момент времени. Мож- но рассмотреть процесс, в котором случайная сила действует на частицу в течение некоторого достаточно продолжительного временного интервала. Для простоты рассмотрим случай, когда в каждый следующий отсчет по времени сила, действовавшая на частицу в предыдущий момент, уменьша- ется в 2 раза. Уравнение для нахождения координаты броуновской частицы (62) не изменится. А для скорости частицы в последовательные моменты времени получим v 1 = 0 + F 0 ∆t m , v 2 = v 1 1 − ∆t mB + F 1 + F 0 2 ∆t m , 47 v 3 = v 2 1 − ∆t mB + F 2 + F 1 2 + F 0 4 ∆t m , v i+1 = v i 1 − ∆t mB + F i + F i−1 2 + . . . + F 0 2 i ∆t m , (74) отсчеты F i , по-прежнему, представляют собой случайные числа. Лабораторная работа 3. Броуновская динамика Целью настоящей работы является исследование динамики броуновской частицы, описываемой уравнением Ланжевена. Содержание работы 1. Выпишите математическую модель для моделирования уравнения Ланжевена. Определите набор входных параметров. 2. Выполните обезразмеривание уравнений и начальных параметров. 3. Выполните программу. 4. Напишите оценку для входных параметров модели, а именно соотно- шение между F max , ∆t и C. Обезразмерьте эту связь и найдите соот- ношение между безразмерными F max и ∆t. 5. Найдите, чему равны средний квадрат скорости и средний квадрат смещения броуновской частицы в случае двумерной модели движения. 6. Сравните полученные результаты с теоретическими. 7. Переделайте программу на случай, когда случайная сила действует продолжительное время (процедура (74)). 8. Напишите выражение для корреляционной функции hF(t)F(t + τ )i в случае, когда случайная сила действует продолжительное время. 9*. Выполните пункты 4–6 для этого случая. 48 Приложение. Комментарии к коду программы Прилагаемая программа на языке С снабжена подробными комментари- ями. Программа реализует моделирование движения броуновской частицы и сбор статистики. Особое внимание следует обратить на функцию langevene_solution, реализующую расчет координат и скоростей для одного шага по времени. Здесь параметр dt задает шаг по времени. В структуре P хранятся коор- динаты и проекции скоростей частицы. Далее следуют функции, которые собирают данные для построения гра- фиков и осуществляют запись и чтение данных о координатах и скоростях в файл. Эти функции можно не изучать подробно. Основной является функция main, где собственно и происходит выпол- нение моделирования и накопление статистики. Переменная count опре- деляет количество проходов. Переменная eps задает точность расчетов. При запуске программы следует ввести число шагов по времени N , шаг по времени ∆t, максимальное значения случайной силы F. Начнется про- цесс моделирования. При нажатии комбинации клавиш Ctrl+C, вычисле- ния приостановятся. Далее следуйте инструкциям на экране. Список литературы [1] Хеерман Д.В. Методы компьютерного эксперимента в теоретической физике. М.: Наука, 1990. [2] Гулд Х., Тобочник Я. Компьютерное моделирование в физике. М.: Мир, 1990. т.1,2. [3] Метод молекулярной динамики в физической химии. Под ред. Ю.К. Товбина. М.: Наука, 1996. [4] Frenkel D., Smit B. Understanding Molecular Simulation. From Algorithms to Applications. San Diego, Academic Press, 2002. [5] Федоренко Р.П. Введение в вычислительную физику. М.: Изд-во Моск. физ.-техн. ин-та, 1994. 49 [6] Allen M.P., Tildesley D.J. Computer Simulation of Liquids, Oxford, Oxford University Press, 1990. [7] Ландау Л.Д.. Лифшиц Е.М. Статистическая физика. Часть 1. М.: На- ука, 1995. [8] Metropolis N., Ulam S. Journal of the American Statistical Association, Vol. 44, No. 247, pp. 335-34 (1949). [9] Metropolis N., Rosenbluth A.W., Rosenbluth M.N., Teller A.H., Teller E., Equation of state calculations for fast computing machines, J. Chem. Phys. 6, 1087 (1953). 50 1 2 |