CC BY
89
УДК 539.3, 519.6
Проблемы применения метода молекулярной динамики при исследовании неравновесных процессов в мезомеханике
И.Ф. Головнев, Е.И. Головнева, В.М. Фомин
Институт теоретической и прикладной механики СО РАН, Новосибирск, 630090, Россия
В настоящей работе проведено исследование динамических и стохастических свойств металлических наноструктур. Показано, что на временных интервалах до 100 пс отсутствует противоречие между законами классической механики и термодинамики, связанное с проблемой обратимости. Доказано, что метод молекулярной динамики, основанный на аналитической механике, может быть использован для обоснования неравновесной термодинамики на микро- и мезоуровнях.
Ключевые слова: мезомеханика, неравновесная термодинамика, динамическая концепция классической механики, метод молекулярной динамики, второе начало термодинамики, стохастические свойства наносистем
Problems of molecular dynamics application in mesomechanical research
on nonequilibrium processes
I.F. Golovnev, E.I. Golovneva, and V.M. Fomin
Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia
The present paper studies dynamic and stochastic properties of metal nanostructures. On time intervals up to 100 ps the laws of classical mechanics and thermodynamics exhibit no contradiction related to a reversibility problem. It is proved that the molecular dynamics method based on analytical mechanics can be used for the substantiation of nonequilibrium thermodynamics on micro- and mesoscales.
Keywords: mesomechanics, nonequilibrium thermodynamics, dynamic concept of classical mechanics, molecular dynamics, second law of thermodynamics, stochastic properties of nanosystems
1. Введение
В серии фундаментальных работ [1-3] на основании анализа экспериментальных и теоретических данных показано, что дальнейшее развитие физической мезоме-ханики и наноматериаловедения должно базироваться на основных положениях неравновесной термодинамики. Но, как известно (см., например, [4]), уравнения неравновесной термодинамики, даже в линейном приближении, могут дать только наиболее общие выводы о неравновесных процессах в твердотельных системах. Более детализированную информацию можно получить, опираясь на соотношения для временных корреляций, которые можно рассчитать для таких систем только численно. Не менее важен и путь прямого численного исследования неравновесных процессов с помощью метода молекулярной динамики. И в первом, и во втором
подходе необходимо проведение динамического расчета, опирающегося на классическую механику. И тут возникает проблема соответствия получаемых численных результатов теоретическим положениям динамической теории как динамической концепции, альтернативной классическому термодинамическому подходу. Наиболее удобным формализмом для анализа этой проблемы является алгебраическая модификация классической механики фон Неймана.
2. Алгебраическая модификация классической механики
Необходимо отметить, что при использовании формализма Лагранжа или Гамильтона динамика механической системы описывается с помощью дифференциальных уравнений. Принципиально иной способ опи-
© Головнев И.Ф., Головнева Е.И., Фомин В.М., 2012
сания эволюции системы базируется на алгебраическом подходе, развитом в работах Купмана [5], фон Неймана [6], Пригожина [7, 8]. (Подробное изложение этого подхода можно найти, например, в работах [9-11].)
В этом подходе рассматривается динамическая система с п степенями свободы с функцией Гамильтона Н в фазовом пространстве Г = {а1, р1,qn, рп} = {у} и множество Р = {/} бесконечно дифференцируемых функций на пространстве Г, называемых наблюдаемыми. При этом уравнения Гамильтона порождают коммутативную группу преобразований пространства Г в себя: ; Г ^ Г! 0{^ = С^Сг., которые, в свою оче-
редь, порождают семейство преобразований алгебры наблюдаемых в себя: Ц; Р ^ Р! йtf (у) = А(у) = f (0{у).
Аналогично квантовой теории вводится понятие состояния механической системы, которое состоит в том, что при многократном повторении опыта получаем вероятностное распределение для всех наблюдаемых. Не делая различия между аналитической механикой и статистической физикой, вводится единое понятие математического ожидания наблюдаемой в состоянии |а):
{А “) = / f (^) Ра (У)^>
т.е. состояние |а) в механике считается заданным, если задано вероятностное распределение ра. При этом вводится понятие чистого состояния:
N
Ра(У) = И8(4 - 4- ^))§(Р — Р- (t))
-=1
е шаоаша! ёаё ё^аШ шйд1у1еу й дашдаааёа1ёа1, 1оёё^а^йё1йу 10 Шё5аааа1ёу 8-бо1ё6ёё.
Эволюция механической системы рассматривается как зависимость от времени среднего значения наблюдаемой в состоянии |а). При этом используются два способа такого описания. Первый — так называемая картина Гамильтона. Вся зависимость от времени перенесена на саму наблюдаемую:
А = {н А} дра = о
И-~{н •А }'~дГ
тогда А | а) = I \и1А(4о^ Ро)]Ра (4о^ Ро)й4оФо.
Другой способ (картина Лиувилля) получается в результате замены переменных: 0(у0 ^ у. Тогда (А | а{) = = IА(4, Р)Рt(4, Р)^4йр. Следовательно, в картине Лиувилля вся зависимость от времени перенесена на вероятностное распределение ра:
А = о дра = {н Р }
Эквивалентность картин А | = ^ А | а^ следует из са-
мого вывода эволюции среднего значения для наблюдаемой.
Нетрудно видеть, что традиционно в аналитической механике используются картина Гамильтона и чистые состояния, а в статистической физике рассматриваются смешанные состояния и картина Лиувилля.
В рамках гамильтоновой механики уравнение движения для наблюдаемой f имеет вид:
А = —{н, А} = -ьт dt
при условии, что f не зависит явно от времени. Необходимо отметить, что именно этот случай наиболее важен на практике. Так, в качестве f могут выступать, например, обобщенные координаты, импульсы, энергия системы и т.д.
В этом выражении {Н, /} — скобки Пуассона, а Ь — оператор Лиувилля, введенный впервые И. Пригожи-ным [7, 8]:
Ь = у А—.дН А
г=1 4 ЭР- дР- %
В этих работах оператор был взят с мнимой единицей для подчеркивания аналогии с уравнением Шредин-гера. Однако в наших дальнейших приложениях это не понадобится и мнимая единица опущена (см., например, [9]).
В случае обобщенно-консервативных систем, когда
дН = о, дt
т.е. связи являются склерономными, а внешние поля — стационарными, общее решение уравнения движения для наблюдаемой можно представить в виде:
А ^) = А (а- ^X р1 ^)) = ъМ^Ао = йц А
Здесь /о = А(а-(о), р- (о)), а 0^) — классический про-пагатор, оператор Грина [9], или оператор эволюции. Известно, что оператор 1)^) имеет вид:
^ (— 1)п#п
) = УЬп,
=о п!
(1)
где Ь = {...{Н,...},...} — п раз вложенные скобки Пуассона.
Таким образом, этот формализм, рассматривающий с единых позиций задачи механики и статистической физики, отвечает исходной постановке задачи: расчет неравновесных процессов в твердых телах методами классической механики. При этом основой описания динамики системы являются операторы эволюции Ц.
Дифференциальные уравнения Гамильтона в рамках формализма фон Неймана переходят в операторные соотношения эволюционного вида:
4(0 = U(t)q(0)• р(0 = и(!) р(о). (2)
При этом все физические положения, заложенные в формализме Гамильтона, сохраняются и в алгебраическом подходе фон Неймана.
) В цитируемой литературе показано, что множество Ц образует однопараметрическую непрерывную группу Ли, т.е. выполняются следующие условия: 1) существует единичный элемент и (о) = 1; 2) существует обратный элемент и_1(0 = и(—т.к. и(—)Ц^) = ); 3) про-
изведение двух элементов также принадлежит этому множеству и(^ )и(^) = и(^ )Ц(^ = и(^ + ^) (ассоциативность умножения).
Эти свойства определяют и динамическую концепцию классической механики.
1. Обращение времени равносильно движению по той же траектории в исходное состояние, т.е. при инверсии времени ) оператор эволюции переходит в обратный: 1,0 (0 = и (-) = и_1(0.
2. Уравнения механики инвариантны не относительно обращения времени (инверсии времени ) : t ^ -), а относительно обращения времени-импульсов ) р ^ ^ — р ^ —р).
Это отчетливо видно для систем с функцией Гамильтона
Н = 1^ + и (41,..., дп) + УЩ (а- ),
-=12ш,- г=1
которая наиболее часто встречается в приложениях. Здесь и(а1,..., 4п) — энергия взаимодействия атомов между собой, а (а- ) — взаимодействие атомов с внешними потенциальными полями. Уравнения Гамильтона, описывающие динамику этой системы, имеют вид:
■ _ Р±_ ■ _ ди(д1; дп) ЪШ(д,)
Ч.Ї , рі ^ ^ .
т Щі дд
(3)
Если провести операцию обращения времени Iг (г ^ -г), то уравнения примут вид:
&, _-Рі,Р, _№(дд' ■'Чп> + ^■2 (4)
дд д%
и не совпадают с уравнениями Гамильтона (3), т.е. они не инвариантны относительно обращения времени. В то же время, если провести обращение времени-импульса 1{ (г ^ -г, р ^ -р), то уравнения остаются инва-
риантными. Подробный обзор работ, посвященных данной теме, можно найти в монографии [12].
При инверсии 1{ оператор Лиувилля приобретает знак «-»: Ь ^ -Ь, а оператор эволюции
ехр{-Ь} ^ ехр{-(-г)(-Ь)} _ ехр{-Ь} остается инвариантным. При исследовании динамики системы удобно провести операцию инверсии I(р в некоторый момент г0. Тогда получаем:
■і _ и(Ч)Ь (г0) _ и(г0)и(г0)д,(0) _
_ ии(2г0)д,(0) _ ■ (2а
Рі _ и(г0)(-Рі(г0)) _ -и(г0)и(г0)Р,(0) _
_-(/(2г0) Р,(0) _-Р, (2г0).
Отсюда следует, что координаты в момент г _ 2г0 остались теми же, а у импульсов остались инвариантными модули, сами же они приобрели противоположный знак. Следствием инвариантности эволюции координат относительно I, „ является то, что и динамика концент-
I, р
рации частиц п(г, (), и динамика полной потенциальной энергии и (х1(, ),..., xN ^)) инвариантны относительно ) . Так как кинетическая энергия системы зависит от квадратов модулей импульсов, то и кинетическая энергия, и локальная температура также инвариантны относительно ) . Следовательно, перенос энергии и массы в макросистеме инвариантны относительно обращения времени-импульса.
Нетрудно видеть, что обобщенные уравнения переноса, описывающие неравновесные процессы [13]: дп + дп(ук) = о
дг дхк
дп{иі) дп{иі) _ дЖ дг дхк дх, ’
д 81 . гк 9 д 81 . 2 .9 дЖ, . -п{V2) | + —| -п<V V) |_ п—<ук),
дхк
(5)
дг I 2
дхк і 2
также инвариантны относительно операции ) р (,г = -, , р-=—иг-, г-= г), как и уравнения классической механики. Действительно, эти уравнения получены из цепочки Боголюбова, при условии парности взаимодействий. Следовательно, они и должны иметь ту же симметрию, что и уравнение Лиувилля для полной системы. На практике используются уравнения, в которых обобщенный тензор напряжения аппроксимируется феноменологическим приближением:
р<Де,.Аик) = Щк -Пк, где Р — статическое давление, а
п ,к _п 'ік +п "к,
п,к _а
8
ди, дик 2 дщ
3 дх1
9
8,,
(6)
дТ
(7)
• +
дхк дх
# ЪТЛ] ~
П к =ПдХ 8-к
дх1
и теплового потока:
Ьк=1 р<д^ 2д^к )=-к^.
2 дхк
Инвариантность уравнений (5) относительно инверсии времени-импульса приводит к аналогичной инвариантности выражений (6) и (7) и, следовательно, определяет симметрию коэффициентов вязкости а, п и теплопроводности к по отношению к оператору ) .
Таким образом, к численным результатам, полученным в рамках метода молекулярной динамики, с помощью которых проводится исследование неравновесных термодинамических процессов, необходимо предъявлять три дополнительных требования: 1) при обращении времени-импульсов в начальный момент времени фазовые траектории исходной системы и системы, подвергнутой инверсии, совпадают; 2) при обращении времени ) система должна вернуться в начальную точку по старой траектории; 3) второе начало термоди-
намики должно выполняться при обращении времени-импульса, т.е. симметрия неравновесной термодинамики совпадает с симметрией уравнений механики.
3. Численная схема
Для анализа поставленной проблемы логично было бы использовать численную схему, полученную в рамках подхода фон Неймана. В качестве примера ниже рассмотрим получение численной схемы для обобщенно-консервативных систем из общего вида оператора эволюции. Численный метод базируется на выражении
f (г + т) _Ц/(т) f (г),
или
~ (-1)птп
/ (г+Т) _Е ^— Ьп/(г).
п_0 п!
Используя свойства оператора эволюции, непрерывное отображение (1) заменяем дискретным отображением
f (г) _ [й(^)Гf (0),
где t = тт. Для построения численной схемы оператор эволюции £/(т) для достаточно малых временных шагов т аппроксимируется оператором Ш1 (т):
Подставляя это выражение в (11) и далее в (10), получаем:
I (-1)птп „
и (т) _Е ьп
(8)
п=о п!
Верхний предел суммирования определяется как величиной т, так и необходимой точностью расчета. Можно было бы использовать в расчетах непосредственно выражение (8), но логичней показать, что наиболее используемая в расчетах методом молекулярной динамики скоростная модификация схемы Верле [14] является следствием (8) и использовать именно известную схему.
Действительно, первое уравнение для координат является просто операторным разложением (2) до второго порядка:
х(к+1) _ хк + — Рк +— Рк т 2т
(9)
Второе уравнение схемы получается из оператор-
(10)
ного разложения для импульсов
р (к+1) = р (к)-1 Ьр(к) + — Ь2 р(к >
1! 2!
и действия операторов Лиувилля на импульсы:
и к)=№;(к>=- Р (к),
V дх ; ’ (11)
Ь2р(к > = Ь( Ьр(к >) = Ь(-р (к >).
С другой стороны, операторное разложение для силы до первого порядка имеет вид:
Р (к+1) = Р (к) - тЬР(к
Следовательно,
~ Р(к+1) - Р(к )
- ЬР(к) = Р-------Р—.
Р( к+1) _ Р(к)+1 (р (к+1)+р (к)),
(12)
т.е. второе уравнение скоростной модификации Верле.
Везде ниже в численных результатах использовался целый ряд обезразмеривающих множителей: координаты измерялись в 1о-1о м, время — в 1о-13 с, масса — в 1о-27 кг. Исходя из этих значений, получены следующие единицы измерения: для энергии — в 1о-21 Дж, скорости — 1о3 м/с, силы — 1о-11 Н, давления — 1о9 Па.
4. Физическая система
Численную иллюстрацию поставленной задачи наиболее уместно провести для металлических систем. Во-первых, это один из основных объектов исследования мезомеханики. Во-вторых, в настоящее время практически отсутствуют систематизированные работы по изучению стохастических свойств многочастичных металлических потенциалов, аналогичные работам, посвященным моделям потенциала твердых сфер или лен-нард-джонсовских систем. Так, для систем твердых сфер основные закономерности были получены еще группой Б. Олдера при создании самого метода молекулярной динамики [15] и в дальнейшей серии работ (см., например, [16-18]). Наиболее полные исследования динамических и стохастических свойств систем твердых сфер проведены в серии работ группы В.Я. Рудяка [19, 20].
Свойствам систем частиц с потенциалами взаимодействия типа Леннард-Джонса уделено еще больше внимания, и списки работ достигают нескольких сотен. Приведем здесь лишь одни из самых первых работ этого направления [21-25].
Работы [26, 27] посвящены исключительно исследованию стохастических свойств систем с потенциальным взаимодействием Леннард-Джонса. Подробный обзор можно найти, например, в работах [28-32]. В то же время исследованию динамических и стохастических свойств твердотельных структур посвящены лишь немногие работы. Так, влиянию внешних шумов на процессы в системе посвящены работы [8, 9, 33, 34]. Реализация стохастических граничных условий при молекулярно-динамическом моделировании была предложена в работе [35]. Исследование влияния граничных воздействий на стохастические характеристики одномерных цепочек проведено с помощью метода молекулярной динамики в работе [36]. В работах [37, 38] на основании обобщения статистики Больцмана-Гиббса исследована термодинамика структурно-скейлинговых переходов при пластической деформации твердых тел.
Таким образом, проблема соответствия получаемых численных результатов с помощью метода молекуляр-
т
ной динамики и основных динамических концепций для твердотельных (металлических) систем оказалась открытой. Для поставленной задачи не столь важна универсальность потенциала, сколько то, что он описывает основные свойства многочастичных металлических систем и имеет аналитический вид, что позволяет избежать потери точности при интерполировании потенциала и сил по табличным численным значениям.
В настоящее время значительный прогресс в численных расчетах свойств металлов обусловлен появлением потенциалов, рассчитанных из первых принципов, например, на основе метода внедренного атома (МВА) [39-45]. Разработан и хорошо себя зарекомендовал в описании тепловых и механических свойств переходных металлов так называемый «квантовый потенциал Саттона-Чена» [46]. Кузнецов В.М., Каминский П.П. и др. (см., например, [47]) развили метод модельного функционала электронной плотности для получения многочастичных потенциалов межатомного взаимодействия в металлах.
Однако в связи с поставленной задачей использовался аналитический МВА-потенциал для меди, предложенный в [48, 49]. Это позволило, с одной стороны, изучать свойства металлических систем, а с другой — не терять точность при интерполировании функций по точкам, как например, при использовании потенциала [50].
Все основные численные иллюстрации будут представлены для следующих систем.
1. Кластеры меди сферической формы с начальным расположением атомов в узлах кристаллической ГЦК-ячейки идеального кристалла. Начальные импульсы полагались равными нулю. В связи с тем, что в кластерах значительную роль играет поверхность, данная система относится к неравновесным [51-55].
2. Эти же кластеры, но «охлажденные» с помощью метода искусственной вязкости [56] до криогенных температур, а затем нагретые до нужных температур с помощью метода стохастических сил [52]. Эти системы находятся в термодинамическом равновесном состоянии.
3. Кластер, находящийся в начальный момент времени в термодинамическом равновесном состоянии, с заданной начальной скоростью сталкивается с одномерным потенциальным барьером. При этом барьер моделируется отталкивательной ветвью потенциала Лен-нард-Джонса
и =еь V] (х > 0).
В связи с этим начальная скорость кластера отрицательна. В результате неупругого столкновения кластер переходит в термодинамическое неравновесное состояние и улетает от потенциального барьера.
4. Кластер в форме прямоугольного параллелепипеда разогрет неравномерно, т.е. приведен в неравновесное состояние, в результате чего можно проследить эволюцию прихода в термодинамическое равновесное состояние в результате процесса теплопередачи.
Радиусы кластеров менялись от 1 до 4 нм, а число атомов — от 369 до 22 663 соответственно, т.е. отличие составляло почти два порядка. Это позволило провести исследование зависимости динамических и стохастических характеристик от числа атомов в достаточно широком интервале размеров наноструктур.
Необходимо также отметить, что все представленные системы относятся к замкнутым и изолированным, т.е. они не обмениваются массой и энергией с внешней средой. Это позволило проверить наличие стохастич-ности металлических структур, т.е. зарождение хаоса, обусловленное свойствами самой системы.
5. Закон сохранения энергии
Очевидно, что без иллюстрации выполнения закона сохранения энергии в численных результатах и правильной работы численной схемы (п. 3) не может быть и речи об исследовании более тонких свойств системы. В качестве примера на рис. 1 приведена зависимость логарифма ошибки по энергии
8
Ег - Ег 0
г0
9
•100%
от времени для трех значений временного шага т = 0.01, 0.001 и 0.0001 для эволюции кластера с радиусом 1 нм и температурой Т = 50 К.
Как видно, уменьшение шага по времени на порядок приводит к уменьшению ошибки на два порядка, что доказывает, что эта схема является схемой второго порядка точности. Кроме того, используя шаг т = 0.001 (10-16 с), мы имеем систему с постоянной энергией, флуктуирующей около постоянного значения с ошибкой порядка 10-8 %. Это позволяет сделать вывод о том, что численное моделирование удовлетворяет физическому условию замкнутости системы.
Поскольку ряд динамических характеристик зависит от числа атомов системы, дополнительно проводилось исследование зависимости ошибки по энергии от этой величины. При этом было показано, что при увеличении числа атомов на два порядка ошибка не изменилась.
Также исследовалось влияние температуры на ошибку по энергии. Показано, что она не изменяется в исследуемом интервале температур и составляет примерно 10-8 % для шага по времени т = 0.001 (10-16 с).
Не менее важен для исследуемой проблемы и ответ на вопрос: Долго ли результаты численной схемы будут соответствовать условию изолированности системы? На рис. 2 приведена зависимость ошибки на интер-
вале до 100 пс (на два порядка больше предыдущего примера). Как видно, энергия остается постоянной с той же точностью. Кроме того, необходимо отметить, что времена релаксации в твердотельных наноструктурах значительно меньше этого интервала.
На этом же рисунке (рис. 2) приведено влияние не-равновесности процесса на величину ошибки энергии. За эталон берется также криогенное состояние (линия 1). Линия 2 соответствует ошибке при отражении кластера, находящегося в начальный момент времени при криогенной температуре, от потенциального барьера с начальной скоростью 10 м/с. Как видно, ошибка не изменилась, хотя столкновение с барьером сопровождается неравновесными процессами в кластере.
В то же время не столь значительное отклонение от равновесия (линия 3), когда в начальный момент половина кристалла охлаждена до криогенной температуры, а половина — нагрета до 50 К, привело к росту ошибки на четыре порядка. Еще более значительный рост ошибки, почти на шесть порядков, обусловил выбор в качестве начального состояния координаты и импульсы первой модели (линия 4), когда атомы кластера располагались в узлах идеальной кристаллической ячейки, а импульсы полагались равными нулю.
Этот результат показывает, что на ошибку сильно влияет отклонение от термодинамического равновесия начального состояния системы. В то же время численные результаты, полученные с помощью метода молекулярной динамики, отвечают физической постановке задачи — рассмотрению эволюции изолированной системы.
6. Инверсия времени и времени-импульса в методе молекулярной динамики
Проверка инвариантности уравнений классической механики относительно инверсии 1г р импульса-времени при использовании численного метода молекулярной динамики проводилась следующим образом. После того как были получены начальные данные для координат и импульсов {г- , р -} для перечисленных выше вариантов задач, готовились начальные данные для второй траектории {г1- , р1-}, причем г1- = г , а у импульсов изме-
нялся знак: р1- =-р I. Расчет двух траекторий проводился параллельно в одной и той же программе (для уменьшения машинных округлений чисел при запоминании массивов данных), но для обращенной траектории шаг по времени изменялся на противоположный: т1 = -т. В отличие от (9), (12), расчетная схема приобрела вид для координат:
х(к+1 = хк - — Рк + — Рк
-РК +-
т 2т
и для импульсов:
р(к+1) = р(к) -1 (Р(к+1) + р(к)).
В процессе расчета сравнивались как сами траектории {г| , р|}, так и все основные механические характеристики, являющиеся функциями координат и импульсов системы. В процессе расчета было показано, что разница между этими характеристиками равна в точности нулю для эволюции всех процессов, перечисленных выше, в связи с чем нет возможности привести их графическую иллюстрацию. Следовательно, численные результаты, полученные в рамках метода молекулярной динамики для металлических систем, не противоречат основному динамическому принципу классической механики — инвариантности динамических уравнений относительно инверсии времени-импульсов.
Вторая операция инверсии времени в приложениях вызывает больший интерес, т.к. позволяет провести исследование времени обратимости [20], которое в этой работе определяется как максимальное значение времени движения по фазовой траектории, из которого при инверсии времени 1г можно вернуться в заданную 8-окрестность начальной точки. Близкое по физическому смыслу имеет понятие времени обратимости тобр, введенное в работах Ю.А. Кравцова [57, 58]. Он отождествлял время обратимости с временем замкнутости системы, в течение которого на систему не действуют внешние флуктуационные поля. Близкий физический смысл имеет и время «горизонта предсказуемости», введенное в работе [59] как интервал времени, за который малое возмущение начальных условий для динамической переменной или же малое изменение параметров системы приводит к непредсказуемости ее дви-
Рис. 1. Зависимость ошибки от времени для значений шага по времени т = 0.0001 (1), 0.001 (2), 0.01 (3)
Рис. 2. Зависимость ошибки от времени для различных неравновесных процессов в системе (т = 0.001, R = 1 нм)
жения. Также близкий смысл имеет и «время динамической памяти» т характеризующее время забывания начальных условий, введенное в работах [60-62].
Учитывая все особенности этих определений относительно задачи, поставленной в настоящей работе, — возможности использования приближенных численных решений для обоснования теоретических положений, в работе были проведены численные эксперименты для перечисленных выше процессов следующим образом. После формирования начальных данных для всех систем, описанных выше, проводилось численное интегрирование «вперед» и после определенного интервала времени (составляющего, как правило, для твердотельных наноструктур 10-12 с) проводилось запоминание массивов координат и импульсов всех атомов (в бесформатном режиме). На втором этапе решения задачи проводилось интегрирование «назад»: в цикле считывались из памяти массивы координат и импульсов, «оставленные» при интегрировании «вперед», и интегрирование шло с отрицательным шагом по времени (инверсия времени 1{). В процессе расчета выводилась зависимость механических характеристик от времени старта «назад» при движении к начальной точке. При этом проводилось и сравнение отклонений этих величин в начальной точке траектории.
В качестве примера на рис. 3 приведено отличие полной энергии исходной и «пришедшей» назад системы в начальной точке от длины интервала возвращения в зависимости от температуры равновесного кластера (кривые 1, 2) и для неравновесного начального состояния (модель 1, кривая 3). Видно, что неравновесность приводит к снижению точности примерно на порядок, но ошибка все равно составляет 10-11 %. Кроме того, видно, что эта точность не зависит от времени возврата, т.е. система не забывает о начальном состоянии в течение 10-10 с.
Расчеты зависимости времени памяти от размеров (для температуры кластера 100 К) были проведены для кластеров с размерами 1, 1.5 и 2 нм, что соответствует изменению числа атомов примерно на порядок. При этом было показано, что размеры системы на обратимость или время памяти не влияют.
Аналогичные расчеты были проведены для столкновения кластера с потенциальным барьером и показано, что точность совпадения в начальной точке также порядка 10-12 а система имеет время памяти большее 50 пс. Для этого процесса также имеется высокая обратимость по времени, а с увеличением начальной температуры в кластере при столкновении с барьером наблюдается даже уменьшение ошибки по энергии.
Условие, что время обратимости по механическим характеристикам для наноструктур более 100 пс, привело к тому, что для них оказалось невозможным реализовать метод, предложенный в работе [19], — усред-
Рис. 3. Зависимость ошибки в полной энергии от времени обращения для различных состояний системы (т = 0.001, R = 1 нм): Т = 25 (1) и 100 К (2), неравновесный кластер (3)
нение по временным интервалам, большим времени обратимости как усреднение для представителей ансамбля. По этой причине для характеристики отклонения траекторий от начальных точек при возвращении в начальную точку ниже рассмотрен ряд характеристик, аналогичных корреляциям, но отличающийся от них тем, что они не усреднялись по ансамблю или по времени.
Наиболее прозрачный смысл имеют 8-окрестности начальной точки, определяемые как (для сравнения с другими авторами ниже введена окрестность в подпространстве скоростей)
= 11 Л2
8г ^узк^1 Г|0
11 V, л2
8В =—ЛтттЕ(V -vi0) .
(13)
(14)
Эти величины обезразмерены на характерные параметры задачи: Д — радиус кластера, а = ЗкТ/т —
среднеквадратичная максвелловская скорость (аналогично работе [26]). Здесь г|0 и Уг0 — координаты и скорости атомов в первоначальный момент времени, а гг- и уг- — координаты и скорости обращенных траекторий в начальной точке.
В качестве примера на рис. 4 приведены значения ^ 8Ъ. (Для величины 8Г получены те же по величине результаты.) Для неравновесной модели (линия 1) температура при определении юТ взята равной 50 К, полученной в качестве средней величины в процессе эволюции к термодинамически равновесному состоянию. Линии 2 и 3 соответствуют обращению времени для эволюции кластеров с температурой 25 и 100 К.
Как видно, наблюдается слабый экспоненциальный рост 8-окрестностей как для равновесных, так и для неравновесных состояний, хотя это сказывается только на временных интервалах до 100 пс, а само отклонение при этом меньше 10-5 %. Это позволяет сделать вывод о наличии эффекта забывания начальных данных в системе, хотя и пренебрежимо малом на временных интервалах, представляющих интерес для кинетики наноструктур. Полученное время обратимости для твердотельных наноструктур примерно в 100-500 раз больше
аналогичной характеристики для аргона в работах [19, 20]. По-видимому, при значительном увеличении интервала интегрирования время обратимости будет проявляться и для механических характеристик системы.
Необходимо отметить, что аналогичные исследования столкновения с барьером показали, что угол наклона в зависимости ег-окрестности сильно уменьшился, а е-окрестность вообще перестала зависеть от времени обращения.
Таким образом, эволюция рассмотренных выше систем, рассчитанная в рамках численного метода молекулярной динамики, на рассмотренном интервале времени в 100 пс с большой точностью подчиняется классическому динамическому принципу возврата в начальное состояние при обращении времени.
7. Локальная неустойчивость в металлических наноструктурах
Локальная устойчивость в настоящей работе рассмотрена в соответствии с теорией Ляпунова аналогично методике, примененной в работах [19, 20]. Численные исследования проводились следующим образом. В начальный момент времени имелись начальные координаты г и скорости Vг- всех атомов системы, найденные в качестве начальных данных, как было описано выше. Далее, с помощью датчика случайных чисел находятся координаты и импульсы возмущенной системы в 8-окрестности начального положения системы в соответствии с выражениями
х(1) = х- + £8еУ-1 = •••>
= V х + ^т,
где £ — случайное число в интервале (-1; 1), а 8е — расстояние между траекториями в начальный момент. Далее интегрирование для двух траекторий идет параллельно в одной программе и через определенное число шагов по времени находится расстояние ег и £в между траекториями в координатном и импульсном подпространствах по формулам (13), (14).
В качестве примера на рис. 5 приведена зависимость логарифма расстояния ег от времени (для начального расстояния ^ ег = -2.23) для Т = 25 К. Видно, что с увеличением размера наноструктур, зависимость расстояния между фазовыми траекториями от времени ослабевает и для кластеров с радиусом более 2 нм практически отсутствует, т.е. по формальным признакам траектории локально устойчивы, и лишь для самых малых размеров имеется очень слабая локальная неустойчивость, что не совпадает с результатами для газов в работах [19, 20, 26, 27]. Правда, в работах [26, 27] локальная неустойчивость исследовалась для одной и той же системы, но две траектории получались при использовании разных расчетных шагов в численной схеме, что иллюстрирует, скорее, свойство самой схемы. Однако, хотя
-14 -)------;-------1------1-------1-------
0 20 40 60 80 1. ПС
Рис. 4. Зависимость е^ от времени обращения для разных начальных данных системы (т = 0.001): кластер в неравновесном состоянии (1), Т = 25 (2) и 100 К (3)
в настоящей работе локальная неустойчивость исследовалась по методике, предложенной в работах [19, 20], результаты показывают отличие свойств газовых и твердотельных систем.
Общим для всех вариантов расчетов настоящей работы, как видно из рис. 5, является экспоненциальный рост расстояния между фазовыми траекториями для времен меньших 1 пс. Особенно это явление ярко видно для расстояний в импульсном подпространстве (рис. 6) на малых временах релаксации. Это обусловило проведение детализированного исследования процессов в системах для начальных интервалов времени.
В качестве примера на рис. 7 приведена зависимость температуры в основной и дополнительной системе от времени. Видно, что в возмущенной системе протекает неравновесный процесс, а температура в течение 1 пс релаксирует к некоторой постоянной величине около 140 К.
Для анализа наличия термодинамического равновесия в системе использовался критерий Колмогорова [63] для вероятности реализации гипотезы распределения Максвелла по скоростям для атомов в системе. При этом было показано, что система релаксировала к термодинамическому равновесному состоянию на интервале времени до 0.5 пс. На этом же интервале наблюдается и наиболее значительный рост расстояния между фазовыми траекториями систем.
1дег
■3.0 -)-----------------1---------------1---------------1---------------1---------------
0 20 40 60 80 1, ПС
Рис. 5. Зависимость расстояния между фазовыми траекториями в координатном подпространстве от времени. Температура Т = 25 К, радиус кластера R = 1 (1), 1.5 (2), 2 (3), 3.5 нм (4)
0.00 0.02 0.04 ^ пс
Рис. 6. Зависимость расстояния между фазовыми траекториями в импульсном подпространстве от времени. Т = 100 К, R = 1 нм
Таким образом, численные исследования показали, что при приготовлении начальных данных для дополнительной системы в е-окрестности основной системы с помощью метода Монте-Карло формируется система в неравновесном термодинамическом состоянии, а полная энергия этой системы отличается от исходной энергии. Экспоненциальный рост расстояния между фазовыми траекториями систем соответствует релаксации дополнительной системы в равновесное состояние, а время роста является временем релаксации. В то же время после релаксации системы в равновесное состояние расстояние между траекториями остается практически постоянным.
Данный вывод привел к необходимости проведения численного эксперимента расчета расстояния между траекториями, исходно находящимися в термодинамическом равновесном состоянии, причем их начальные данные даже могут и не находиться в малой е-окрестности. Для иллюстрации этого положения на рис. 8 приведена зависимость расстояния еу для двух термодинамически равновесных систем при температуре 25 и 50 К. (Для величины ег результаты аналогичны.) Как видно, расстояние между траекториями остается постоянным. Аналогичные результаты были получены в расчетах систем при одинаковых температурах 25 К, но в разных микросостояниях.
Отсюда следует, что фазовые траектории, рассчитываемые в рамках молекулярно-динамической модели для металлических наноструктур, удовлетворяют дина-
Т к ,_____________________________________
о -3------------------------1----------------------1----------------------1---------------------1---------------------
0 0 0.2 0.4 0 6 0 8 I ПС
Рис. 7. Зависимость температур в основной (1) и в дополнительной системе в е-окрестности (2) от времени. Т = 25 К, R = 1 нм
Рис. 8. Зависимость расстояния между фазовыми траекториями в импульсном подпространстве от времени, Т = 50 К
мическому критерию локальной устойчивости по Ляпунову на временном интервале до 100 пс.
8. Перемешивание фазовых траекторий для металлических наноструктур
Свойство перемешивания выражается как свойство расцепления временны х корреляций некоторых динамических переменных f и g
Нш R(/г, g) = 0.
Т
В работе рассмотрена автокорреляционная функция скоростей
R (г) = 11ш 1 / Е (v - (г + £), V- (£)). (15)
Т^~ г 0 N -=1
Ниже приведена иллюстрация для случая кластера сферической формы с радиусом 2 нм (число атомов в системе N = 2 899). В качестве примера на рис. 9 приведена зависимость функции
<v(г), v(0)> =1Е (V-(I), V-(0))
N -=1
от времени.
Видно, что уже на временах порядка 3 пс функция осциллирует около нуля. Исследование зависимости < v(t + £), V(г)> от сдвига по времени £ в разные моменты времени t показало, что на временах до 0.25 пс они совпадают с большой точностью, а слабое отличие наблюдается в районе «хвостов», где они также осциллируют около нуля и вклад в интеграл (15) не дают. На рис. 10 приведена зависимость относительной величины Я(г)/<V(0)V(0)> от времени. Видно, что корреляция на интервале порядка 1 пс становится равной нулю, т.е. система обладает свойством перемешивания или условием конечности времени релаксации. Необходимо подчеркнуть, что примерно то же самое значение временного интервала было получено выше (см. п. 7) при исследовании релаксации системы в термодинамическое равновесное состояние.
9. Проблема Ферми-Паста-Улама для металлических наносистем
Проблема перемешивания траекторий, рассмотренная выше, близка к проблеме стохастизации или тепловой диссипации механического возмущения кристалла, возникшего под влиянием внешней среды. Впервые эта
(у®, v(0))
0.008 0.000
-0.008
0 1 2 ^ пс
Рис. 9. Зависимость < v(^), v(0)> от времени. Радиус кластера 2 нм, Т = 25 К
проблема была поднята Дебаем в связи с исследованием теплопроводности в твердых телах, а позже была сформулирована Э. Ферми, С. Уламом и Дж. Пастой [64]. Они попытались дать объяснение перераспределению энергии между колебательными модами и термализации системы осцилляторов, связанных нелинейными силами. Более подробно эта проблема рассмотрена в работах [65, 66]. В связи с этим используемый метод молекулярной динамики должен быть проверен на задаче термализации.
В качестве примера была рассмотрена задача соударения кластера с потенциальным барьером. Начальная температура кластера была близка к нулю для снятия тепловых флуктуаций, а начальная скорость — 10 м/с. Кластер налетает на барьер, находящийся в начале координат, в связи с этим у него начальная скорость отрицательна, т.к. направление движения в сторону уменьшения координат.
На рис. 11, а приведена зависимость компоненты скорости центра масс Усх от времени. Видно, что соударение начинается около момента времени 10 пс (скорость начинает заметно изменяться), а прекращается около 30 пс. На рис. 11, б представлена зависимость изменения внутренней потенциальной энергии кластера и хаотической (тепловой) доли кинетической энергии атомов от времени. После прекращения взаимодействия кластера со стенкой продолжается энергообмен в колебательных степенях свободы кристалла еще примерно
ВД)/<у(0), v(0))
0.4 0.2 0.0 -0.2
0.0 0.2 0.4 0.6 0.8 I, пс
Рис. 10. Зависимость относительной величины корреляции от времени. Радиус кластера 2 нм, Т = 25 К
в течение 10 пс. Этот процесс сопровождается тепловыми осцилляциями и постепенным выравниванием потенциальной и хаотической составляющей внутренней энергии, что является признаком установления термодинамического равновесного состояния.
Таким образом, используемый метод молекулярной динамики и многочастичный потенциал Джонсона описывают известный экспериментальный факт термализации.
10. Метод молекулярной динамики и второе начало термодинамики
После того как подвергнуты анализу численные результаты, получаемые для металлических нанокластеров, представляет интерес проверить прямым численным исследованием выполнение второго начала термодинамики. Для этого наиболее удобно проанализировать формулировку Клаузиуса, согласно которой теплота не может переходить сама собой от более холодного тела к более горячему.
В качестве численной иллюстрации ниже приводятся данные расчета процесса теплопередачи в кластере Си, имеющем форму прямоугольного параллелепипеда (число кристаллических ячеек — 20, 3, 3 вдоль соответствующих осей координат) и с начальным распределением температуры, имеющим вид ступеньки. Средние значения начальных температур в левой половине кристалла Т10 = 49.3 К, а в правой — Т20 = 0.04 К. Расчет всех характеристик проводился до установления термо-
\/сх, км/с
ДЕ, 10~14 эрг
Рис. 11. Зависимость х-компоненты скорости центра масс от времени (а) и зависимость логарифма изменения внутренней потенциальной энергии кластера (линия 1) и хаотической доли кинетической энергии атомов (линия 2) от времени (б)
Рис. 12. Зависимость температуры подсистем от времени
динамического равновесия. На рис. 12 приведена зависимость средних значений температур в подсистемах от времени.
Линия 1-2-3 — зависимость температуры от времени в более холодной части кластера, а линия 4 - 53 — в более горячей части кластера. Момент времени 0 соответствует началу процесса теплообмена. Как видно, температуры, в соответствии со вторым началом термодинамики, выравниваются, а время релаксации составляет примерно 20 пс. После этого проводился второй численный эксперимент. Конечные координаты и импульсы атомов (точка 3) брались в качестве начальных данных, и далее решалась система уравнений (4), полученная операцией инверсии из уравнений Гамильтона. При этом делалось столько же шагов «назад». Температуры при этом эволюционировали по своим траекториям также «назад», и система вернулась в первоначальное состояние с температурным профилем в виде «ступеньки». При этом отклонения в данных по траекториям отдельных атомов не превышали 0.01 %, и кривые для обратной эволюции совпали с линиями 1-2-3 и 4-5-3.
Третий численный эксперимент проводился для анализа влияния обращения времени-импульса, при котором уравнения Гамильтона остаются инвариантными, на описание процесса теплопередачи. Система эволюционировала до 4 пс (рис. 12, точки 2 и 5), в соответствии с уравнениями Гамильтона (3), а затем делалось обращение времени-импульсов ,р и система эволюционировала дополнительно 16 пс. Кривые 1-2 и 4-5 претерпели излом (время «пошло назад»), и в системе теплообмен продолжал совершаться совершенно симметрично «положительному» направлению времени (кривые 2- 6 и 5 - 6 соответственно). Видно, что обращение времени-импульса не только не изменяет уравнений движения, но и правильно описывает неравновесный процесс теплопередачи.
Количественной мерой процесса необратимости процессов является энтропия. В работе она рассчитывалась на основе термодинамического определения количества тепла, подведенного или отведенного от под-
систем 8QІ = сіёТ-, и с использованием приближения квазистатического процесса и связи изменения энтропии с изменением количества тепла в подсистемах SQІ• = Т . В результате полное изменение энтропии может быть выражено через изменение температур и пропорционально величине
М = М11п
8 т 9 8 т2 9
1 + ЖЛп 2
Т о 2 <Т20 =
(16)
Здесь с — теплоемкости подсистем; N — количество атомов в подсистемах; Т0, Т — начальные и текущие температуры подсистем.
Постановка численного эксперимента описана выше при анализе теплопередачи. Параллельно расчету температур проводилось вычисление приведенной энтропии (рис. 13). Путь 1-2-3 соответствует прямому решению уравнений движения атомов, путь 3 - 2-1 — инверсии времени, а путь 2-4 — обращению времени импульсов 1г р в момент 4 пс. Как и в случае с температурой, при инверсии времени 11. система вернулась в исходное состояние по этой же линии 3-2-1 (с точностью 0.01 %), т.е. энтропия, как и все другие характеристики, ведет себя обратимым образом по отношению к инверсии времени. В то же время из симметрии графиков 2-3 и 2-4 следует инвариантность энтропии относительно операции инверсии времени-импульса
Кр •
Таким образом, как уравнения механики, так и уравнения, описывающие динамику неравновесных процессов, инвариантны относительно обращения времени-импульсов I р, а эволюция системы и на микро-, и на макроуровне является необратимой. Следовательно, между механическим и кинетическим способом описания системы нет противоречия, а результаты, полученные с применением метода молекулярной динамики в рамках формализма классической механики, позволяют правильно описать необратимые процессы на мезо- и макроуровнях.
11. Заключение
Таким образом, показано, что численные результаты, полученные с помощью метода молекулярной дина-
Рис. 13. Зависимость энтропии от времени
мики для металлических наноструктур, удовлетворяют основным динамическим концепциям классической механики:
- обращение времени равносильно движению по той же траектории в исходное состояние;
- инвариантность уравнений механики не относительно обращения времени, а относительно обращения времени-импульсов;
- выполнение второго начала термодинамики при обращении времени-импульса, т.е. наличие симметрии неравновесной термодинамики относительно инверсии времени-импульса и совпадение с симметрией уравнений механики;
- метод молекулярной динамики на интервалах до 100 пс с большой точностью описывает замкнутые и изолированные наноструктуры.
Это позволило проверить ряд стохастических свойств металлических наноструктур:
- время обратимости для твердотельных наноструктур превышает 100 пс, т.е. система не забывает о начальном состоянии на этом временном интервале;
- фазовые траектории, рассчитываемые в рамках молекулярно-динамической модели, удовлетворяют динамическому критерию локальной устойчивости по Ляпунову на временном интервале до 100 пс;
- временные корреляции на интервале порядка 1 пс становятся равными нулю, т.е. система обладает свойством перемешивания или условием конечности времени релаксации.
Все это позволяет сделать вывод о возможности использования численных результатов, полученных в рамках метода молекулярной динамики, для исследования неравновесных термодинамических процессов и дальнейшего обоснования физической мезомеханики и наноматериаловедения.
Отдельные разделы настоящей работы выполнены в рамках междисциплинарных интеграционных проектов СО РАН №№ 72 и 78.
Литература
1. Панин В.Е., Егорушкин В.Е. Физическая мезомеханика и неравно-
весная термодинамика как методологическая основа наноматериаловедения // Физ. мезомех. - 2009. - Т. 12. - № 4. - С. 7-26.
2. Панин В.Е., Егорушкин В.Е., Панин А.В. Эффект каналирования пластических сдвигов и нелинейные волны локализованной пластической деформации и разрушения // Физ. мезомех. - 2010. -Т. 13. - № 5. - С. 7-26.
3. Панин В.Е., Егорушкин В.Е. Деформируемое твердое тело как нелинейная иерархически организованная система // Физ. мезомех. -2011. - Т. 14. - № 3. - С. 7-26.
4. Терлецкий Я.П. Статистическая физика. - М.: Высшая школа,
1994.- 350 с.
5. Koopman B.O. Hamiltonian systems and transformations in Hilbert space // Proc. Natl. Acad. Sci. U.S. - 1931. - V. 17. - P. 315-318.
6. Neumann J.V Allgemeine Eigenwerttheorie Hermitischer Funktional-
operatoren // Math. Annalen. - 1929. - V. 102. - P. 49-131.
7. Brout R., Prigogine I. Statistical mechanics of irreversible processes. Part VII: General theory of weakly coupled systems // Physica. -1956.- V. 22. - P. 621-636.
8. Balescu R., Prigogine I. Irreversible processes in gases. II: The equations of evolution // Physica. - 1959. - V. 25. - P. 302-323.
9. Балеску P. Равновесная и неравновесная статистическая механика. - М.: Мир, 1978. - Т. 1. - 405 с.
10. Zwanzig R.W. Statistical mechanics of irreversibility // Lectures in Theoretical Physics. V III. - N.Y.: Interscience Publishers, Inc., 1961.
11. Фаддеев Л.Д., Якубовский О.А. Лекции по квантовой механике. -Л.: Изд-во ЛГУ, 1980. - 200 с.
12. Хайтун С.Д. Механика и необратимость. - М.: Янус, 1996. - 448 с.
13. Уленбек Дж., Форд Дж. Лекции по статистической механике. -М.: Мир, 1965. - 307 с.
14. Swope W.C., Andersen H.C., Berens PH., Wilwson K.R. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters // J. Chem. Phys. - 1983. - V. 76. - P. 637-649.
15. Alder B.J., Wainwright T.E. Studies in molecular dynamics. I. General method // J. Chem. Phys. - 1959. - V. 31. - No. 2. - P. 459-466.
16. Alder B.J., Wainwright T.E. Velocity autocorrelations for hard spheres // Phys. Rev. Lett. - 1967. - V. 18. - P. 988-990.
17. Alder B.J., Wainwright T.E. Studies in molecular dynamics. III. Transport coefficients for a hard-spheres fluid // J. Chem. Phys. - 1970. -V. 53. - No. 10. - P. 3813-3826.
18. Alder B.J., Wainwright T.E. Decay of the velocity autocorrelation function // Phys. Rev. A. - 1970. - V. 1. - P. 18-21.
19. Рудяк В.Я., Иванов Д.А. Компьютерное моделирование конечного числа взаимодействующих частиц // ДАН ВШ. - 2003. - Т. 1. -С. 30-38.
20. Рудяк В.Я., Иванов Д.А. Динамические и стохастические свойства открытой системы конечного числа упруго взаимодействующих частиц // Труды НГАСУ. - 2004. - Т. 7. - № 30. - С. 47-58.
21. Rahman A. Correlations in the motion of atoms in liquid argon // Phys. Rev. - 1964. - V. 136. - No. 2. - P. 405-411.
22. Verlet L. Computer «experiments» on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules // Phys. Rev. - 1967. -V. 159. - No. 1. - P. 98-103.
23. Лагарьков А.Н., Сергеев В.М. Вычисление коэффициентов переноса плотных газов и жидкостей методом молекулярной динамики // ТВТ. - 1970. - Т. 8. - С. 1309-1311.
24. Лагарьков А.Н., Сергеев В.М. Исследование переносных и термодинамических свойств аргона методом молекулярной динамики // ТВТ. - 1973. - Т. 11. - С. 513-522.
25. Levesque D., Verlet L., Kurkijarvi J. Computer “experiments” on classical fluids. IV. Transport properties and time-correlation functions of the Lennard-Jones liquid near its triple point // Phys. Rev.
A. - 1973. - V. 7. - P. 1690-1700.
26. Норманн Г.Э., Стегайлов В.В. Стохастические свойства молекулярно-динамической ленард-джонсовской системы в равновесном и неравновесном состояниях // ЖЭТФ. - 2001. - Т. 119. - № 5. -С. 1011-1020.
27. Norman G.E., Stegailov VV. Stochastic and dynamic properties of molecular dynamics systems: Simple liquids, plasma and electrolytes, polymers // Comp. Phys. Comm. - 2002. - V. 147. - No. 1-2. - P. 678683.
28. Лагарьков А.Н., Сергеев В.М. Метод молекулярной динамики в статистической физике // УФН. - 1978. - Т. 125. - № 3. - С. 409448.
29. Билер Дж. Машинное моделирование при исследовании материалов. - М.: Мир, 1974. - 416 c.
30. Хеерман Д.В. Методы компьютерного эксперимента в теоретической физике. - М.: Наука, 1990. - 176 c.
31. Alien M.P., Tildesley DJ. Computer Simulation of Liquids. - N.Y.-Oxford: Clarendon Press, 1987. - 385 p.
32. Метод молекулярной динамики в физической химии. - М.: Наука, 1996. - 334 с.
33. Хорстхемке В., Лефевр Р. Индуцированные шумом фазовые переходы. - М.: Мир, 1987. - 400 с.
34. Лихтенберг А., Либерман М. Регулярная и стохастическая механика. - М.: Мир, 1984. - 528 с.
35. Dmitriev A.I., Psakhie S.G. Stochastical Boundary Conditions for Computer Simulation of Materials under Loading of Different Type // Proc. II China-Russia Symposium on Advanced Materials and Processes, Xian, China. - 1994. - P. 152-157.
36. Сараев Д.Ю. Исследование нелинейного отклика в твердом теле на атомном уровне при высокоскоростном нагружении: Дис. ... канд. физ.-мат. наук. - Томск, 1998. - 129 c.
37. Наймарк О.Б., Баяндин Ю.В., Леонтьев В.А., Пермяков С.Л. О термодинамике структурно-скейлинговых переходов при пластической деформации твердых тел // Физ. мезомех. - 2005. - Т. 8. -№ 5. - С. 23-29.
38. Пантелеев И.А., Frostey C., Наймарк О.Б. Структурно-скейлин-говые переходы и универсальность статистики флуктуаций при пластическом течении металлов // Вычислительная механика сплошных сред. -2009. - Т. 2. - № 3. - С. 70-81.
39. Foiles S.M. Calculation of the surface segregation of Ni-Cu alloys with the use of the embedded-atom method // Phys. Rev. B. - 1985. -V. 32. - No. 12. - P. 7685-7693.
40. Foiles S.M., Baskes M.I., Daw M.S. Embedded-atom method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys // Phys. Rev. B. - 1986. - V 33. - P. 7983-7991.
41. Jacobsen K.W., Norskov J.K., Puska M.J. Interatomic interactions in the effective-medium theory // Phys. Rev. B. - 1987. - V 35. - P. 74237442.
42. Nelson J.S., Sowa E.C., Daw M.S. Calculation of phonons on the Cu (100) surface by the embedded-atom method // Phys. Rev. Let. -1988.- V. 61. - No. 17. - P. 1977-1980.
43. Daw M.S., Baskes M.I. Embedded atom method derivation and application to impurities, surfaces, and other defects in metals // Phys. Rev.
B. - 1984. - V. 29. - No. 12. - P. 6443-6453.
44. Daw M.S., Foiles S.M. Summary abstract: Calculations of the energetic and structure of Pt (110) using the embedded atom method // J. Vac. Sci. Technol. A. - 1986. - V. 4. - No. 3. - P. 1412-1413.
45. Daw M.S. Calculation of the energetic and structure of Pt (110) reconstruction using the embedded atom method // Surf. Sci. - 1986. -V. 166.- No. 2-3. - P. L161-169.
46. Cagin T., Dereli G., Uludogan M., Tomak M. Thermal and mechanical properties of some fcc transition metals // Phys. Rev. B. - 1999. -V. 59. - P. 3468-3473.
47. Физическая мезомеханика и компьютерное конструирование материалов: В 2 т. / Под ред. В.Е. Панина. - Новосибирск: Наука,
1995. - Т. 1. - 298 с.; - Т. 2. - 320 с.
48. Johnson R.A. Analytic nearest-neighbor model for fcc metals // Phys. Rev. B. - 1988. - V. 37. - P. 3924-3931.
49. Johnson R.A. Alloy models with the embedded-atom method // Phys. Rev. B. - 1989. - V. 39. - P. 12554-12559.
50. Voter A.F., Chen S.P. Accurate interatomic potentials for Ni, Al, and NisAl // Mat. Res. Soc. Symp. Proc. - 1987. - V. 82. - P. 175-180.
51. Болеста А.В., Головнев И.Ф., Фомин В.М. Исследование процесса соударения сферического кластера меди с жесткой стенкой методом молекулярной динамики // Физ. мезомех. - 2000. - Т. 3. -№ 5. - С. 39-46.
52. Болеста А.В., ГоловневИ.Ф., Фомин В.М. Плавление на контакте при соударении кластера никеля с жесткой стенкой // Физ. мезо-мех. - 2001. - Т. 4. - № 1. - С. 5-10.
53. ГоловневИ.Ф., Головнева Е.И., Фомин В.М. Молекулярно-динамическое исследование столкновения нанокластеров друг с другом и с подложкой // Физ. мезомех. - 2007. - Т. 10. - № 2. - С. 5-13.
54. Головнев И.Ф., Головнева Е.И., Фомин В.М. Расчет термодинамических свойств наноструктур методом молекулярной динамики // Физ. мезомех. - 2007. - Т. 10. - № 5. - С. 71-76.
55. ГоловневИ.Ф., Головнева Е.И., Фомин В.М. Молекулярно-динамическое исследование поверхностного натяжения наноструктур // МТТ. - 2010. - № 3. - С. 45-55.
56. Головнева Е.И., Головнев И.Ф., Фомин В.М. Моделирование ква-зистатических процессов в кристаллах методом молекулярной динамики // Физ. мезомех. - 2003. - Т. 6. - № 6. - С. 5-10.
57. Кравцов Ю.А. Фактические границы гипотезы замкнутости и классические парадоксы кинетической теории // ЖЭТФ. - 1989. -Т. 96. - № 5(11). - С. 1661-1665.
58. Кравцов Ю.А. Случайность, детермированность, предсказуемость // УФН. - Т. 158. - № 1. - С. 93-122.
59. Lighthill J. The recently recognized failure of predictability in Newtonian dynamics // Proc. Roy. Soc. Ser. A. - 1986. - V. 407. - P. 35-50.
60. Норман ГЭ. Стохастизирующий фон молекулярной динамики. Уравнения молекулярного движения. Необратимость // Сб. научн. сообщ. V Всесоюз. конф. по строению и свойствам металлических и шлаковых растворов. Ч. I. Теория жидких и аморфных металлов. - Свердловск: УНЦ АН СССР, 1983. - С. 58-62.
61. Заславский Г.М. Стохастичность динамических систем. - М.: Наука, 1984. - 272 c.
62. Заславский ГМ. Гамильтонов хаос и фрактальная динамика. -М.-Ижевск: НИЦ «Регулярная и хаотическая динамика», Ижевский институт компьютерных исследований, 2010. - 472 с.
63. КолмогоровА.Н. Теория вероятностей и математическая статистика. - М.: Наука, 1986. - 534 с.
64. Ферми Э. Научные труды. Т. 2. Сер. «Классики науки». - М.: Наука, 1972. - 717 с.
65. Заславский Р.М., Сагдеев Р.З. Введение в нелинейную физику. От маятника до турбулентности и хаоса. - М.: Наука, 1988. - 368 с.
66. Тода М. Теория нелинейных решеток. - М.: Мир, 1984. - 264 с.
Поступила в редакцию 18.05.2012 г., после переработки 14.09.2012 г.
Сведения об авторах
Головнев Игорь Федорович, к.ф.-м.н., снс ИТПМ СО РАН, golovnev@itam.nsc.ru Головнева Елена Игоревна, к.ф.-м.н., снс ИТПМ СО РАН, elena@itam.nsc.ru Фомин Василий Михайлович, академик, директор ИТПМ СО РАН, fomin@itam.nsc.ru
