Извлечение колебательного спектра из расчёта BOMD

[Vit Nhoc]

Я хочу реализовать в Chemcraft такую опцию: открыть файл с расчётом BOMD (отдельной молекулы), выбрать “BOMD tools/Show vibrational spectrum” и чтобы программа построила ИК или КР спектр. Наверно, пользователю нужно будет ещё указать точку, где начинается эргодичность.
Подскажите, как подступиться к этому делу:
1) Молекулярная механике в BOMD отдельной молекулы – она скорее ньютоновская, или скорее шредингеровская, или на каких-то костылях?
Если ньютоновская – может вообще из этой идеи ничего не выйдет? Если не ошибаюсь, когда считают колебательный спектр обычными методами (по силовому полю), решают вековое уравнение и это вековое уравнение подразумевает именно квантовую, а не классическую механику (ядерную).
2) Как вообще происходит поглощение молекулой кванта ИК-излучения? Какими факторами определяется интенсивность полосы в ИК спектра?
3) Если в BOMD молекула “гуляет” для фиксированной полной энергии (кинетической плюс потенциальной), не приведёт ли предлагаемый подход к тому, что в построенном спектре будут только малые частоты, а для больших частот молекула не будет “допрыгивать” (т.е. кинетическая энергия молекулы в BOMD намного меньше собственной энергии моды).
4) Как совместить то, что колебательный спектр квантуется, с тем, что в BOMD ничего на уровне ядер не квантуется, ядра просто гуляют по ППЭ. Этот вопрос скорее дополняет первый.

[madschumacher]

Наверно, пользователю нужно будет ещё указать точку, где начинается эргодичность.

Само собой!

1) Молекулярная механике в BOMD отдельной молекулы – она скорее ньютоновская, или скорее шредингеровская, или на каких-то костылях?

Молеkулярная механиkа -- это отдельная песня. Под ней подразумеваются (последние лет 20) параметризованные эмпиричесkие силовые поля (CHARMM, AMBER и иже с ними). А kонkретно BOMD -- это решение kлассичесkих уравнений (чаще всего Ньютона) на B-O ППЭ (Kэп ). Это если эмулировать NVE ансамбль. Если же нужны свойства NVT ансамбля, то добавляются еще термостаты (их kуча всяkих, это отдельная тема). Но в любом случае это будет kлассичесkая траеkтория (NVE или NVT) в фазовом пространстве.
Делается очень просто -- даются начальные условия (чаще всего генерятся случайно из Маkсвелловсkого распределения), и потом интегрируют 2й заkон Ньютона F=-grad V=ma. Т.е. на kаждом шаге траеkтории считается градиент ППЭ. И... собственно всё. Вот kлассичесkая MD траеkтория и получена!
Справедливости ради, надо отметить, что есть определенные термостаты, kоторые позволяют получить kвантовые NVT распределения (например, т.н. Colored Noise Thermostates, реализованы в CP2K и может где ещё).

Если ньютоновская – может вообще из этой идеи ничего не выйдет?

Выйдет!

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

Именно таk, но разницы ниkаkой -- в гармоничесkом приближении kвантовая и kлассичесkая задачи имеют одно и то же решение (не в смысле траеkтории или волновые фунkции), а спеkтр, т.е. собственные числа и собственные веkтора (а.k.а. частоты и формы нормальных kолебаний). Kcт, интенсивности тоже вроде при этом совпадают (насчет этого не очень уверен ).

2) Как вообще происходит поглощение молекулой кванта ИК-излучения? Какими факторами определяется интенсивность полосы в ИК спектра?

Ооой, ну там вроде множитель | <n| d |m> |² стоит в первом (дипольном) приближении (если в kвантовом случае рассматривать задачу поглощения, для перехода n -> m, или наоборот ). А в kлассичесkом (вроде? ) через флуктуационно-диссипационную теорему (+ теорему Винера-Хинчина) выписывается основное расчетное выражение (т.е. спеkтр = Фурье-образ автоkорреляционной фунkции неkоторой величины, в случае ИK -- дипольного момента). Ну это и логично, т.k. этот спеkтр = "слепоk" взаимодействия системы с внешним элеkтромагнитным полем (a.k.a слабым внешним возмущением).

3) Если в BOMD молекула “гуляет” для фиксированной полной энергии (кинетической плюс потенциальной), не приведёт ли предлагаемый подход к тому, что в построенном спектре будут только малые частоты, а для больших частот молекула не будет “допрыгивать” (т.е. кинетическая энергия молекулы в BOMD намного меньше собственной энергии моды).

В целом это здравое соображение. И это действительно таk (если под собственной энергией моды имелась в виду ZPE). В лучшем случае -- там действительно у высоkих частот существенно ниже интенсивность, но это решают при помощи немалого зоопарkа kвантовых мультиплиkативных фунkций--поправоk для спеkтра. Самая популярная (вроде) -- это ~ ν * tanh (hν/2kT). И есть еще kуча (со всяkими разными свойствами). Но бывает это не помогает из-за т.н. "flying ice cube effect" (если термостат хреновый ) -- в этом случае из высоkочастотных мод энергия чересчур сильно сливается в низkочастотные. Тут уж, kонечно, ничего не поделать уже. Но по частям его (спеkтр) все равно можно смотреть.
Там загвоздkа вообще есть другая и еще хуже -- ангармонизм для высоkих частот сильно недооценивается (т.е. линии ближе k гармоничесkим, чем k реальным). Но с этим, вроде можно бороться повысив температуру и вообще это уже эффеkт несkольkо другого порядkа малости по частоте.

4) Как совместить то, что колебательный спектр квантуется, с тем, что в BOMD ничего на уровне ядер не квантуется, ядра просто гуляют по ППЭ. Этот вопрос скорее дополняет первый.

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

P.S. BOMD не единственный вариант ab initio MD. Есть еще методы типа CPMD, и иже с ними, где элеkтронная задача сопряжена (неадиабатичесkи) с ядерной.
P.P.S. А для чисто kвантовых свойств есть PIMD + техниkи для уменьшения необходимого числа реплиk (CN термостаты для того и лепили изначально), а там еще Троттеровckое разложение, и мне недавно поkазали статью, что можно с multiple time steps тоже это сделать. И там то спеkтры уже совпадают по определению (ну с точностью до всяkих численно--реализационных эффеkтов).

[Vit Nhoc]

Молеkулярная механиkа -- это отдельная песня.

Сорри, оговорился - я имел в виду молекулярную динамику.

Именно таk, но разницы ниkаkой -- в гармоничесkом приближении kвантовая и kлассичесkая задачи имеют одно и то же решение

А в ангармоническом варианте - разница будет большая?

Не очень понятно, что такое термостат и зачем он нужен, если в BOMD и так задаётся какая-то полная энергия системы, соответствующая очевидно какой-то энергии. Или термостат - это для конденсированной фазы?

Я пока не понял принцип, как можно попробовать посчитать ИК-спектр. Предположим, у нас посчитано 1000 шагов в BOMD. Может быть, надо взять 1000*1000 комбинаций, для каждой посчитать энергию перехода (это разность между потенциальными энергиями на данных шагах) и вероятность перехода (это как-то через дипольные моменты на этих шагах). Эти 1000*1000 комбинаций в сумме дадут полный спектр с уширением. Или я всё не так понимаю?

[madschumacher]

А в ангармоническом варианте - разница будет большая?

Для высокочастотных колебаний (по моему скромному опыту, ~ выше 2000 см^-1) до порядка 100 см^-1 в худшем случае. Но люди как-то живут с этим...

Не очень понятно, что такое термостат и зачем он нужен, если в BOMD и так задаётся какая-то полная энергия системы, соответствующая очевидно какой-то энергии. Или термостат - это для конденсированной фазы?

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

Я пока не понял принцип, как можно попробовать посчитать ИК-спектр. Предположим, у нас посчитано 1000 шагов в BOMD. Может быть, надо взять 1000*1000 комбинаций, для каждой посчитать энергию перехода (это разность между потенциальными энергиями на данных шагах) и вероятность перехода (это как-то через дипольные моменты на этих шагах). Эти 1000*1000 комбинаций в сумме дадут полный спектр с уширением. Или я всё не так понимаю?

Совершенно не так. Сначала нужно посчитать автокорреляционную функцию дипольного момента <d(t), d(0)> : <d(t_i), d(0)> = Σ_j (d(t_j) * d(t_j+i) ), (ну там вроде еще нормировка нужна). А потом делается Фурье-преобразование полученной зависимости (для скорости лучше использовать FFT). Это будет "сырой" спектр f(ν). Потом этот спектр f(ν) умножается на корректирующую функцию q(ν) (например, q(ν) ~ ν * tanh (hν/(2kT)) ), которая правит интенсивности. В результате финальный спектр будет I(ν) = f(ν) * q(ν).
Там еще есть некоторые аспекты численной реализации, типа того, что преобразование Фурье для конечной функции дает кривую обрыва (по электронографическому прошлому может Вы это еще помните ), поэтому тут, как и в электронографии для построения кривой радиального распределения, используются оконные функции (например, та же электронографическая exp(-a*s²) -> exp(-(t/τ)²) = Гауссовское окно, τ - полуширина окна).
И для NVE динамики вообще рекомендация использовать не одну траекторию, а много, и с разными начальными условиями (т.к. каждая из них вообще не эргодическая ни в одном глазу), а потом полученные спектры усреднить.

P.S. в автокорреляционной функции я мог немного напортачить, но суть именно такая.

[Vit Nhoc]

Поясните, как строить автокорреляционную функцию.
У нас есть функция dx(t), dy(t), dz(t) - компоненты дипольного момента молекулы от времени t, которое начинается с 0 или с to (эргодичность).

Строим функцию:



Рисунок взял из википедии, здесь f(t) - это dx(t), dy...

Интегрировать придётся не от минус бесконечности до плюс бесконечности, а от t0 до tmax

Мы получим три функции - для трёх компонентов дипольного момента. Что с ними делать дальше?

Ещё мне непонятно: в интеграле стоит f(t-T); если T например находится в середине BOMD-траектории, то эта f будет иметь отрицательные и положительные аргументы, а у нас же нет дипольного момента для отрицательного времени?

У вас в вашей формуле вместо f(t-T) стоит d(j+i), тогда отрицательного аргумента не будет, зато будет аргумент больше tmax.

[madschumacher]

По-моему что-то вот такого типа должно быть:
Если я нигде не ошибся. По-сути это сумма трех автокорр. функций (для dx, dy и dz). Отрицательных времен нет -- и не надо.

[madschumacher]

Ну вот что дальше...

[Vit Nhoc]

Прошу прощения за навязчивость.
1) <d(t),d(0)> - зачем записывать функцию от времени в такой сложной форме, что здесь делает d(0)
2)tanh- это тангенс?
3) в середине второй картинки есть интеграл по времени, в который входит <d(t),d(0)>. Куда здесь пропало то что дипольный момент - это не одно число, а вектор? Ведь на выходе получаем I(v) - это число.

[madschumacher]

1) <d(t),d(0)> - зачем записывать функцию от времени в такой сложной форме, что здесь делает d(0)

Обозначение -- общепринятое. (см, например,

DOI:

http://dx.doi.org/10.1021/acs.jpcb.5b03323

).
Я от него тоже не очень в восторге, откровенно говоря.

2)tanh- это тангенс?

Гиперболический тангенс (в русскоязычной литературе обозначается часто как th(x) )
https://ru.wikipedia.org/wiki/%D0%93%D0 ... 0%B8%D0%B8

3) в середине второй картинки есть интеграл по времени, в который входит <d(t),d(0)>. Куда здесь пропало то что дипольный момент - это не одно число, а вектор? Ведь на выходе получаем I(v) - это число.

Автокорреляционная функция дипольного момента -- это скаляр (сумма автокорреляционных функций x, y и z компонент). Соотв. <d(t),d(0)> пишется только для краткости записи. на самом деле:
<d(t),d(0)> = <d_x(t),d_x(0)> + <d_y(t),d_y(0)> + <d_z(t),d_z(0)>

[Vit Nhoc]

Я пробую запускать расчёты BOMD. Возник вопрос, как поднять колебательную энергию молекулы. Я запустил BOMD ацетона и атомы только болтаются туда-сюда, но нет вращения групп CH3 (т.е. не преодолевается барьер на ППЭ вращения).

[madschumacher]

А какая прога?
Какой вариант MD -- NVE или NVT?
Обычно там в качестве начальных условий генерируются скорости из Максвелловского распределения случайно => можно:
а) запустить новую траекторию (если там с seed-ом все в порядке),
б) сгенерировать начальные скорости в какой-нибудь нормальной MD проге (типа TINKER), или уж самому написать программку (вообще это и в Excel можно сделать ручками ).
+: Вы достаточно траектории нагенерировали, чтобы сказать, что торс. колебаний нет? Они же на ~ 2 порядка медленнее bond stretch колебаний!
++: А температура симуляции вообще достаточна для преодоления барьера? Эксп. барьер (http://www.sciencedirect.com/science/ar ... 526590144X) ≈365 K => при 300 K классически в NVE варианте (если в этом приближении барьер такой же или выше эксп.) Вы его заведомо точно не преодолеете.

[Vit Nhoc]

А какая прога?
Какой вариант MD -- NVE или NVT?
Обычно там в качестве начальных условий генерируются скорости из Максвелловского распределения случайно => можно:
а) запустить новую траекторию (если там с seed-ом все в порядке),
б) сгенерировать начальные скорости в какой-нибудь нормальной MD проге (типа TINKER), или уж самому написать программку (вообще это и в Excel можно сделать ручками ).
+: Вы достаточно траектории нагенерировали, чтобы сказать, что торс. колебаний нет? Они же на ~ 2 порядка медленнее bond stretch колебаний!
++: А температура симуляции вообще достаточна для преодоления барьера? Эксп. барьер (http://www.sciencedirect.com/science/ar ... 526590144X) ≈365 K => при 300 K классически в NVE варианте (если в этом приближении барьер такой же или выше эксп.) Вы его заведомо точно не преодолеете.

Я считаю на Gaussian-е.

Прошу разъяснить разницы между NVE и NVT

> А температура симуляции вообще достаточна для преодоления барьера

Я и хочу узнать, как в расчёте повысить эту температуру симуляции.

Пока я запускаю просто:

#P B3LYP/6-311G(D,P) bomd(maxpoints=4000,ntraj=1)

[madschumacher]

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

nve_system.png

2. закрытые (модель в стат.физике -- канонический ансамбль = NVT)

nvt_system.png

3. открытые (модель в стат.физике -- большой канонический ансамбль = μVT)

mvt_system.png

.
В ab inito MD у нас система -- это молекула (ну или несколько молекул).

md_system.png

Соответственно, если молекула просто ходит в вакууме туды-сюды по своей ППЭ с постоянной изначально заданной энергией -- это модель изолированной системы, т.е. NVE ансамбль => NVE MD.

nve_md.png

Чтобы смоделировать закрытую систему -- нам собственно нужна внешняя система и сопряжение по температуре. Для этого в MD используются термостаты -- некоторые эмуляторы внешней среды.

nvt_md.png

Эти термостаты добавляют и забирают энергию у молекулы с помощью обмена энергией (как в закрытой системе прям! ), таким образом моделируя закрытую систему = NVT.

nvt_md_thermostat_action.png

Ура, Товарищи!!

[madschumacher]

Я считаю на Gaussian-е.

Никогда им для этого не пользовался -- к сожалению личного опыта нет. Но есть такой волшееебный сайт http://www.gaussian.com/g_tech/g_ur/k_bomd.htm (называется мануал к Gaussian). И там вроде написано, что температуру можно задать через волшееебный keyword "Temperature"...
Вот только про наличие термостатов в этом пакете я так до сих пор ничего не понял...

Прошу разъяснить разницы между NVE и NVT

См. пред. сообщение.

[Vit Nhoc]

Новые вопросы:
1) В статье которую вы прислали (Fischer, JPC) сначала посчитали 20 коротких траекторий, получили 20 начальных геометрий, и для них просчитали 20 больших траекторий. Зачем это нужно? Может быть, короткие траектории считались с высокими кинетическими энергиями, чтобы молекула свободно ходила по локальным минимумам, а потом уже длинные траектории рассчитывались для комнатной температуры?
2) Если я правильно понимаю, формулы 2) и 3) в статье мне не нужны – это просто расчёт дипольного момента, который Gaussian и так выдаёт.
3) Не очень понятно, почему используется скаляр (автокорреляционная функция), т.е. простая сумма компонент x,y,z, ведь в геометрии длина вектора – это не сумма его компонент x,y,z, а корень из суммы их квадратов.
Нельзя ли ещё ссылок на статьи, где приводятся все эти формулы?

[madschumacher]

1) В статье которую вы прислали (Fischer, JPC) сначала посчитали 20 коротких траекторий, получили 20 начальных геометрий, и для них просчитали 20 больших траекторий. Зачем это нужно? Может быть, короткие траектории считались с высокими кинетическими энергиями, чтобы молекула свободно ходила по локальным минимумам, а потом уже длинные траектории рассчитывались для комнатной температуры?

Я уже не так хорошо помню, но очень похоже на обычную equilibration phase. Там по-моему не особо жуткий конф.состав и лишние издевательства над системой не нужны.

2) Если я правильно понимаю, формулы 2) и 3) в статье мне не нужны – это просто расчёт дипольного момента, который Gaussian и так выдаёт.

Поскольку это статья фактически есть описание того факта, что в nwchem теперь есть bomd, то наверное да (статьи перед глазами нет).

3) Не очень понятно, почему используется скаляр (автокорреляционная функция), т.е. простая сумма компонент x,y,z, ведь в геометрии длина вектора – это не сумма его компонент x,y,z, а корень из суммы их квадратов.
Нельзя ли ещё ссылок на статьи, где приводятся все эти формулы?

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

А ссылки можно искать в самой статье. Если лень - могу послать через неделю.

[ Post made via Android ]

[Vit Nhoc]

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

Извините за настойчивость, но вот посмотрите: если немного повернуть систему координат, то дипольный момент молекулы не изменится, а скалярная сумма его компонентов dx,dy,dz изменится. Или эта автокорреляционная функции должна считаться именно в определённос системе координат?

[madschumacher]

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

Но вращательно - колебательное взаимодействие может вносить немалые вклады, вроде. Но, в норм мд прогах вообще трансляц.и вращения отделяют. Gamess us, прравда этим не страдает. Про гауссиан хщ

[ Post made via Android ]

[madschumacher]

А так вращение само по себе, если оно не отделено, хорошо видно в спектре (ну тоже спектр поглощения же). А трансляции по идее будут давать Доплеровские сдвиги и уширения.

[ Post made via Android ]

[Vit Nhoc]

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

[ Post made via Android ]

Я и говорю, дипольный момент сам по себе (корень из суммы компонент по x,y,z) не меняется при повороте осей.
А сумма этих компонент самих по себе (кстати они берутся со знаками или их модули?) действительно меняется.