Методика моделирования и расчета свойств неравновесных фаз молекулярных систем Текст научной статьи по специальности «Математика»

УДК 519.688+536.4.032.2

doi: 10.18101/2304-5728-2017-1-72-77

МЕТОДИКА МОДЕЛИРОВАНИЯ И РАСЧЕТА СВОЙСТВ НЕРАВНОВЕСНЫХ ФАЗ МОЛЕКУЛЯРНЫХ СИСТЕМ1

© Герман Евгений Иванович

преподаватель кафедры общей физики Бурятский государственный университет Россия, 670000, Улан-Удэ, ул. Смолина, 24а E-mail: net-admin@list.ru

© Цыдыпов Шулун Балдоржиевич

доктор технических наук, заведующий кафедрой общей физики Бурятский государственный университет Россия, 670000, Улан-Удэ, ул. Смолина, 24а E-mail: shulun@bsu.ru

Рассмотрены принципы модернизации метода молекулярной динамики для возможности моделирования неравновесных фаз молекулярных систем. Приведены результаты апробации данной методики в компьютерном эксперименте.

Ключевые слова: молекулярная динамика, фазовые переходы, неравновесные системы, алгоритм Верле.

Введение

Метод молекулярной динамики (МД) является уникальным инструментом моделирования термодинамических систем, позволяющим прослеживать фазовые траектории вещества. В настоящее время существует ряд программных решений основанных на МД, позволяющих моделировать молекулярные системы в равновесных состояниях. Однако, известные программные продукты (CHARMM, LAMMPS, HOOMD, GROMACS, NAMD) более направлены на моделирование состояний с заданными термодинамическими параметрами (ансамблей) и малопригодны для имитации термодинамических процессов протекающих с большой скоростью (быстрое охлаждение, быстрое прессование), поэтому неравновесное моделирование с помощью вышеуказанных программ сводится к инициализации систем в одних только критических областях и точках фазовых переходов. Наибольший же интерес представляет эволюция молекулярной системы в широком температурном интервале в процессе достаточно быстрого охлаждения. Нами предложена методика моделирования процессов быстрого изохорического и изобарического охлаждения молекуляр-

1 Работа выполнена при поддержке РФФИ грант №15-02-08204а и базовой части госзадания МОН РФ грант №3822

ной системы, а также приемы расчета физических свойств материалов по результатам компьютерного эксперимента.

1. Классический метод молекулярной динамики

Методика МД заключается в решении уравнений движения для системы многих частиц. Так, силу, действующую на 7-ю частицу, можно представить в форме суммы векторов

с12г{ 0 ШФ(гь)

т , =-> 1 , (1)

Л ^ <1т

где Ф(г) - потенциал взаимодействия частиц системы, т - масса частицы, г - радиус вектор частицы.

Любой алгоритм, реализующий МД, сводится к наиболее оптимальному по производительности и точности численному интегрированию уравнения (1).Широкое применение имеет алгоритм Верле в скоростной форме, согласно которому координаты и скорости частиц на каждой итерации вычисляются в следующей последовательности:

А/2

?*(/ + Д0 = ?*(0+1;,.(0^ + а,.(0—' (2)

V, (/+—)== V, (0 + - а,. (0Ль (3)

' 2 ' 2

га" аТ

Vi(t + А/) = Vi+ у) + + АОу , (5)

где ги V,. а, - радиус-вектор положения, скорость и ускорение 7-й частицы со отв етств енно.

В отсутствие внешних сил система частиц должна прийти к термодинамическому равновесию с заданными температурой, плотностью и давлением. Проведение эксперимента, моделирующего термодинамические процессы, предполагает модернизацию алгоритма путем введения процедур корректировки скоростей (процесс изменения температуры системы) или размеров ячейки моделирования (процесс изменения плотности системы).

2. Моделирование изобарных процессов

Методика моделирования процесса изохорного охлаждения (нагрева) предполагает довольно простое изменение классического алгоритма путем использования ввода процедуры пропорционального изменения всех скоростей частиц системы через каждые N итераций моделирования.

Более сложным в реализации является проведение моделирования процесса изобарного охлаждения. В этом случае предполагается использование процедур корректировки скоростей частиц и баростатирования. Процедура бароститирования должна выполнять изменение размеров

ячеики моделирования и расстояния между частицами для поддержания постоянного давления системы. На практике применяются два основных метода поддержания давления: баростат Андерсена [1] и баростат Бе-рендсена [2]. Метод Андерсена коренным образом меняет алгоритм Верле путем введения обобщенных координат связанных с размерами системы и скоростью изменения этих размеров, этот алгоритм довольно сложно реализуем, а жесткость баростата не позволяет моделировать процессы быстрого изменения температуры. В работе баростата Берендсена необходимо иметь известное значение изотермической сжимаемости для каждой точки фазовой траектории, что для моделирования неравновесных процессов невозможно.

Рапапорт [3] описывает организацию ансамбля с постоянным числом частиц, давлением и температурой. Здесь используется баростат, корректирующий положения частиц и размеры системы на коэффициент /л, выведенный из лагранжиана (А7'7)-ансамбля:

Р(У)-Р (6)

¡л = 1 +

г.. У ч

где Р - заданное давление, Р(1) - давление, рассчитанное на данной итерации моделирования. Такой коэффициент пригоден для небольших флуктуаций давления и практически не поддерживает постоянным давление в процессе быстрого охлаждения. Однако, мы применили этот метод в циклической процедуре (рис. 1), позволяющей поддерживать давление с точностью 1% даже в случаях с резким перепадом температуры.

Вход в процедуру

Да

Выход из процедуры

Рис. 1. Процедура баростатирования

(7)

Моделирование изобарного охлаждения дает возможность численного определения таких характеристик как изобарная теплоемкость Сри изотермическая сжимаемость /Зт.

Так, получив в результате компьютерного эксперимента температурную зависимость энтальпии Н, можно вычислить изобарическую теплоемкость:

' <Ш_ с!Т у Р

Для определения сжимаемости необходимо получение двух изобар с близкими значениями давления Рг и Р2. По разностям плотностей на фазовых кривых изотермическая сжимаемость определяется следующим выражением:

!3Т= 1 Рп-Рп, (8)

<Р> Р,-Р2

где рп и рР2 - плотности систем под давлением /'/ и Р2 при одинаковых температурах.

3. Апробация

На основе описанной методики разработана программа моделирования молекулярных систем МООХ11[4]. С помощью этой программы нами проведен ряд экспериментов, моделирующих процессы равновесного и неравновесного охлаждения системы частиц аргона.

На рис. 2 представлена температурная зависимость теплоемкости Ср аргона полученная по результатам моделирования равновесного охлаждения аргона. Для сравнения здесь же представлены литературные данные

[5].

Рис. 2. Изобарическая теплоемкость аргона при давлении 5 МПа: о - результаты нашего компьютерного эксперимента, х - литературные данные [5].

Как видно, результаты хорошо согласуются с литературными данными, в областях фазовых переходов при температурах ~55 и -155 К прослеживаются характерные изменения в зависимости СР(Т).

На рис. 3 и 4 представлены результаты расчета сжимаемости по данным моделирования равновесного и неравновесного охлаждения соответственно.

10 60 110 160 210

Рис. 3. Изотермическая сжимаемость аргона при давлении 4 МПа по данным компьютерных экспериментов, моделирующих охлаждение аргона со скоростью 109К/с.

1.2 1

0.8 0.6 0.4 0.2 0

10 60 110 160 210

Рис. 4. Изотермическая сжимаемость аргона при давлении 4 МПа по данным компьютерных экспериментов, моделирующих охлаждение аргона со скоростью 1012К/с.

На графике температурной зависимости сжимаемости (рис. 3), полученной по данным эксперимента с равновесным медленным охлаждением, можно разделить газовую, жидкую и кристаллическую фазы.

Для неравновесного процесса охлаждения со скоростью 1012 К/с температурная зависимость сжимаемости (рис. 4) претерпевает резкое изменение в окрестностях точки ~90 К, возможно, это характеризует переход неравновесный газ (перенасыщенный пар) - неравновесная жидкость (смесь газа с кластерами жидкости). Также можно отметить переход в окрестностях ~30 К, вероятно соответствующий переходу нестабильной жидкости в аморфную фазу.

Заключение

Предложена методика моделирования равновесных и неравновесных молекулярных систем, позволяющая рассчитывать изобарическую сжимаемость и изотермическую теплоемкость. Результаты апробации хорошо согласуются с эмпирическими данными и наглядно демонстрируют фазовые закономерности вещества при равновесном и неравновесном охлаждении.

Литература

1. Andersen Н. Molecular dynamics simulations at constant pressure and/or temperature // The Journal of Chemical Physics. — 1980. — Vol. 72, No. 4. — P. 2384 -2393.

2. Berendsen H. J., Postma J. P. Molecular dynamics with coupling to an external bath // The Journal of Chemical Physics. — 1984. — Vol. 81. No. 8. — P. 3684 -3690.

3. Rapaport D. C. The Art of Molecular Dynamics Simulation / D.C. Rapaport. — Cambridge: Cambridge University Press, 2004. — 564 p.

4. Свид. 2016617783 Российская Федерация. Свидетельство об официальной регистрации программы для ЭВМ. Программа моделирования молекулярных систем MDDX11 / Е.И. Герман; заявитель и правообладатель Е.И. Герман (RU). — №2016615005; заявл. 20.08.2016; опубл. 14.07.2016, Реестр программ для ЭВМ. — 1 с.

5. Thennophysical properties of fluid systems [Электронный ресурс]. — NIST. — Режим доступа: http://webbook.nist.gov/chemistry/fluid/, свободный. — Загл. с экрана.

METHODS FOR MODELING AND CALCULATING THE PROPERTIES OF NONEQUILIBRIUM PHASES OF MOLECULAR SYSTEMS

Evgeniv I German

Lecturer, Department of General Physics

Buryat State University

24a Smolina St., Ulan-Ude 670000, Russia

Shuhm B. Tsydypov

Dr. Sci. (Engineering), Department of General Physics

Buryat State University

24a Smolina St., Ulan-Ude 670000, Russia

The article deals with the principles of modernization of molecular dynamics

method for modeling of nonequilibrium phases of molecular systems. We present

the results of approbating this method in the computer experiment.

Keywords: molecular dynamics, phase transitions, nonequilibrium systems, Verlet's

algorithm.