Отрицательное давление в молекулярной динамике

[Ben]

Всем привет!
Считаю давление в молекулярной динамике через вириал по классической формуле
$PV=Nk_bT+\frac{1}{3}\langle\sum\limits_{i=1}^{N} r_iF_i^{(int)}\rangle$

Для газов все работает (рис. 1), постоянные параметры P,V,T, показаны флуктуации в окрестности точки, которая соответствует уравнению Менделеева-Клаперона. То есть, грубых ошибок в алгоритме расчета давления, вроде, нет.

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

Почему-то мне кажется, что такое поведение давления хорошо известно, подскажите, куда копать?

[realeugene]

Не знаю, что именно вы считаете, но тензор напряжений, например, в растянутой железяке, вполне может иметь отрицательные компоненты на главной диагонали.

[Ben]

realeugene в сообщении #1520097 писал(а):

Не знаю, что именно вы считаете, но тензор напряжений, например, в растянутой железяке, вполне может иметь отрицательные компоненты на главной диагонали.

Формулу же привел... Считаю термодинамическое давление. При увеличении времени усреднения ситуация не меняется. Для жидкости справедлив закон Паскаля, так что пример с железякой в которой присутствуют стационарные "сдавливающие" силы - явно не мой случай.

[Pphantom]

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

[realeugene]

Ben в сообщении #1520109 писал(а):

Для жидкости справедлив закон Паскаля, так что пример с железякой в которой присутствуют стационарные "сдавливающие" силы - явно не мой случай.

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

Я не помню эту модель. Может и не знаком с нею. Но вот сейчас увидел ваши графики - и что-то мне кажется, что у вас вычислительный алгоритм получился неустойчивым, судя по разбросу посчитанных точек. Может и сама модель неадекватна для воды при атмосферногм давлении и комнатной температуре. В лубом случае, жаловаться на отрицательные точки при такой дисперсии странно.

[lel0lel]

растянутая жидкость

[realeugene]

lel0lel в сообщении #1520141 писал(а):

растянутая жидкость

Про максимальные термометры по ссылке - бред. В максимальном термометре (которымм является и медицинский) столбик ртути рвётся, образуется зазор.

[lel0lel]

Может так, ссылок на статьи нет и в англоязычной Вики этой информации нет. Но возможно, что имелось ввиду, что до некоторой температуры всё же не рвётся. Скажем до 35, но в любом случае всё это слова, никаких подтверждений. Ещё как-то читал про ксилемы -- водопроводящие каналы в высоких деревьях. Тоже можно погуглить.
https://en.m.wikipedia.org/wiki/Xylem

[madschumacher]

Ben в сообщении #1520094 писал(а):

Считаю давление в молекулярной динамике через вириал по классической формуле

Во-первых, у Вас формула для МД неправильная. В МД моментальное изотропное давление считается через выражение
$3 P(t)V(t) = 2 ( E_\mathbf{K}(t) - \Xi (t)) $, где $\Xi = -\frac{1}{2}\sum_i \mathbf{F}_i \mathbf{r}_i$ (Ваш вириал), а $E_\mathbf{K}(t) = \sum_i \frac{m_i}{2} v_i^2(t)$ (кинетическая энергия всех частиц). То же выражение, что Вы используете, работает только средних значений уравновешенного газа.
В связи с этим возникают следующие вопросы: как уравновешивалась система, какие термостаты/баростаты использовались, как считался объём (объём периодической ячейки, ячейки с жёсткими стенками, или какое ещё извращение)?

Потом большой вопрос в том, какое у Вас давление и какие у Вас поля? Чем вода моделируется? TIP4P? Квантовая химия? При заданных условиях вода вообще газ?

[sergey zhukov]

realeugene
Да, я совсем недавно первый раз задумался над тем, почему медицинский ртутный термометр всегда показывает последнюю измеренную у человека температуру, даже если он долго лежит в холодном месте. И почему его нужно встряхивать. Если сравнить медицинский термометр с уличным, то можно увидеть, что у медицинского в месте соединения трубки и шарика есть очень сильное сужение трубки. Ртуть проталкивается через него, когда расширяется. А когда начинает сжиматься, то рвется в этом месте из-за сильного трения о стекло. Вряд-ли ртуть в этом термометре находится в растянутом состоянии после охлаждения.

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

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

[Ben]

Цитата:

Во-первых, у Вас формула для МД неправильная....
Потом большой вопрос в том, какое у Вас давление и какие у Вас поля? Чем вода моделируется? TIP4P? Квантовая химия? При заданных условиях вода вообще газ?

Вода моделируется потенциалом TIP3P, молекул 58 штук.
Объем считается динамически, по алгоритму Convex Hull (картинка).
Параметры PVT задавал по разному, последний вариант: термостат: с коррекцией энергии на краях распределения, давление: сила прикладывается к граничным молекулам, выделенным алгоритмом Convex Hull.
Отрицательные давления появляются когда молекулы начинают сбиваться в клубок (то ли капля, то ли кристалл).

Цитата:

Во-первых, у Вас формула для МД неправильная....

Это было-бы самым хорошим вариантом. Кажется, разница у нас только в знаке перед вириалом. На сегодня я разобрался с каждым шагом на пути от теоремы вириала к приведенной формуле, минус перед "внутренним" вириалом возникает, когда переходят от суммы по молекулам к сумме по взаимодействующим парам, а силу записывают через потенциал взаимодействия.

Настольная книга [M.P. Alen, D.J. Tildesley Computer Simulation of Liquids, формула (2.54)]
Eще подробный вывод этой формулы есть у [Kaznesis, Statistical Thermodynamics and Stochastic Kinetics, глава 8.1.1 Application of the virial theorem: equation of state for non-ideal systems].

Я конечно проверю вывод еще раз, но ведь для газовой фазы все работает. Что не так с конденсатом?

Convex Hull:

[DimaM]

Ben в сообщении #1520179 писал(а):

Отрицательные давления появляются когда молекулы начинают сбиваться в клубок (то ли капля, то ли кристалл).

Скорее всего, начинаются колебания объема. Когда капля сжимается, давление подскакивает. Когда расширяется - давление уменьшается и становится отрицательным.
В вашем случае 58 молекул, считай, 2/3 находятся на поверхности, и объемные величины типа давления вообще плохо определены.

[Ben]

Цитата:

Ben в Скорее всего, начинаются колебания объема. Когда капля сжимается, давление подскакивает. Когда расширяется - давление уменьшается и становится отрицательным.

Так и есть! Колебательный режим виден визуально. Но как с этим бороться? Очень хочется построить изотермы в районе фазового перехода.
Колебательный режим дает ровно нулевое среднее давление.
Полно статей, где народ рапортует о том как хорошо совпадает диаграмма состояния для потенциалов TIP3P, TIP4P. Как же они считают давление?

[madschumacher]

Ben в сообщении #1520179 писал(а):

Параметры PVT задавал по разному, последний вариант: термостат: с коррекцией энергии на краях распределения, давление: сила прикладывается к граничным молекулам, выделенным алгоритмом Convex Hull.

А термостат какой? И какие силы прикладываются, в соответствии с каким алгоритмом?

Ben в сообщении #1520179 писал(а):

Кажется, разница у нас только в знаке перед вириалом.

Нет, там нет разницы со знаками, у меня минус в $\Xi$ ещё запихан. Собственно, я возмущался первой частью ($NkT$): Вы туда ставите заданную температуру, или моментальную, высчитанную из кинетической энергии? Записано, как первый вариант, что неверно, но если имелся в виду второй, то всё норм.

Ben в сообщении #1520179 писал(а):

Объем считается динамически, по алгоритму Convex Hull (картинка).

Вообще без периодических условий, строго говоря, объём не определён, поэтому можно позволить себе пляски с бубнами. Например, Вы можете отмасштабировать этот объём каким-нть коэффициентом $s \geq 1$, чтобы учесть, что молекулы воды -- это не точки, а размазанные в пространстве электронные плотности (по-сути добавить аналог solvent accessibility). Тогда флуктуации давления упадут по абсолютной величине.

Ben в сообщении #1520182 писал(а):

Полно статей, где народ рапортует о том как хорошо совпадает диаграмма состояния для потенциалов TIP3P, TIP4P. Как же они считают давление?

Во-первых, Вы можете попробовать нормально определить объём, задав граничные условия. Во-вторых, Вас интересуют не моментальные величины, которые из-за флуктуаций могут быть какими угодно, а средние значения ансамбля, усреднённые по достаточно большому уравновешенному куску траектории. В-третьих, Вы можете всё то же самое проделать с нормальными баростатами, типа Берендсоновского.

[Ben]

Цитата:

Собственно, я возмущался первой частью ($NkT$)

Этот член у меня считается через кинетические энергии, как у Вас написано.

Цитата:

Во-первых, Вы можете попробовать нормально определить объём, задав граничные условия. Во-вторых, Вас интересуют не моментальные величины, которые из-за флуктуаций могут быть какими угодно, а средние значения ансамбля, усреднённые по достаточно большому уравновешенному куску траектории. В-третьих, Вы можете всё то же самое проделать с нормальными баростатами, типа Берендсоновского.

Иерархия усреднений в программе имеется, я уже сам в ней запутался. А вот с баростатом вопрос серьезный, у меня прикладываются силы к граничным молекулам, эти силы могут порождать и поддерживать колебательный режим. Молекулы на границе все время обновляются, получается, я подталкиваю их к центру распределения и жду когда они вернутся обратно.

По быстрому получить результаты - не получилось, спасибо, что подтолкнули мысль в правильном направлении, попробую реализовать что-то из стандартных баростатов.