Расчёт термодинамических функций молекул

[Vit Nhoc]

Я хочу сделать расчёт термодинамических функций (энтропии, энергии Гиббса и пр.) молекул по молекулярным параметрам (структура, колебательные частоты), рассчитанным из квантовой химии. Эти расчёты уже делают программы вроде Gaussian или Orca, но у моей утилиты будет основное применение - когда в молекуле есть очень маленькие частоты, они вносят большую ошибку в энтропию из-за ангармоничности, поэтому их имеет смысл искусственно повысить.
Мне нужны книги, где приведены эти формулы. Есть И. Н. Годнев, "Вычисление термодинамических функций по молекулярным данным". Там я нашёл годные формулы для колебательной энтропии и ZPE. Но с формулами для вращательной энтропии, как мне показалось, есть какая-то путаница: я нашёл формулу для молекулы типа симметричного волчка, у которого есть два момента инерции, но не нашёл формулу для несимметричной молекулы с тремя компонентами Ia, Ib, Ic:



Там далее в книге есть формула для расчёта суммы по состояниям асимметричного волчка, но надо ещё перевести эту сумму по состояниям в энтропию. В связи с этим ещё вопрос - насколько полезно было бы сделать в моей утилите расчёт не только энтропии, ZPE и энергии Гиббса, но и суммы по состояниям? Расскажите, что такое вообще сумма по состояниям и для чего её используют.
Подскажите другие книги с аналогичными формулами. Мне ChatGPT дал такие названия книг, и ничего из этого я не нагуглил:

"Физическая химия" А.И. Китаева и др.;
"Физическая химия: курс лекций" П.В. Балашова и др.;
"Физическая химия. Термодинамика, статистическая физика и кинетика" А.М. Кузнецова и др.

[Vit Nhoc]

ChatGPT дал мне такую формулу для вращательной энтропии линейной молекулы, это бред?

S_rot = R * ln[(T/T0)^(1/2)] * (πI / h^2)^(1/2)

[Гесс]

Вам возможно подойдет
https://wiki.fysik.dtu.dk/ase/ase/therm ... istry.html - это общий обзор
https://wiki.fysik.dtu.dk/ase/_modules/ ... istry.html - это как оно закодено
Да, малые частоты зачастую насильно повышают, но критерий отсечки не определен, люди традиционно используют 8 или 12 обратных сантиметров, но стандарта нет.

[Vit Nhoc]

Вам возможно подойдет
https://wiki.fysik.dtu.dk/ase/ase/therm ... istry.html - это общий обзор
https://wiki.fysik.dtu.dk/ase/_modules/ ... istry.html - это как оно закодено
Да, малые частоты зачастую насильно повышают, но критерий отсечки не определен, люди традиционно используют 8 или 12 обратных сантиметров, но стандарта нет.

Хм, 8 или 12 см-1? А я считал что 100 или 300...

[Гесс]

300 это вообще какой то безумный критерий, вращение метилов в пропане - уже 200.
100 - это дефолт у Shermo, можно читать тут https://s3-eu-west-1.amazonaws.com/item ... ies_v1.pdf
Орка в примерах дает отсечку (CutOffFreq) на 35, интересное обсуждение скажем https://orcaforum.kofo.mpg.de/viewtopic ... req#p40943
Здесь https://aip.scitation.org/doi/10.1063/5.0061187 речь идет о "where ω0 is an user defined rotor cutoff value typically between 20 and 100 cm−1."
Например https://pubs.rsc.org/en/content/article ... d2cy01337a пишет

Since the harmonic approximation often leads to inaccurate entropies for low- frequency modes, all obtained frequencies have been raised to 12 cm−1 if they were below this value, as described in detail in earlier work.42,43

из них ссыль 42 https://pubs.acs.org/doi/full/10.1021/cs5014267

After carefully checking that the associated nuclear motions corresponded largely to translations or rotations of the species, the frequencies were replaced by 12 cm–1. This value is chosen because the entropy contribution from a 12 cm–1 harmonic vibration approximates that from a pseudotranslational/rotational degree of freedom, as discussed in detail in previous work on zeolite-catalyzed alkene–methanol reactions. (23)

Где ссыль 23 https://www.sciencedirect.com/science/a ... via%3Dihub

. For large-amplitude nuclear motions, the harmonic approximation inherently breaks down and numerical noise becomes significant, as expressed by the presence of spurious imaginary values among the frequencies below 100 cm1. After carefully inspecting the related nuclear motions to verify that they corresponded largely to translations or rotations of the species, these frequencies were replaced by 12 cm1. The same was done with real frequencies below 12 cm1. For adsorbates, at most two such frequencies per species were replaced, while for transition states at most one frequency per species was replaced. The value of 12 cm1 is chosen as a crude estimate of the frequency of a harmonic oscillator with an entropy contribution similar to a translational or rotational degree of freedom. Section S2.1 assesses that this approximation is in reasonable agreement with the fully harmonic as well as the mobile adsorbate approaches [46,47]

При этом последняя цитата - это Норсков, это ИМХО авторитет.

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

-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies = 1.000000000 (already applied!)

0: 0.00 cm**-1
1: 0.00 cm**-1
2: 0.00 cm**-1
3: 0.00 cm**-1
4: 0.00 cm**-1
5: 0.00 cm**-1
6: 0.64 cm**-1
7: 0.72 cm**-1
8: 0.72 cm**-1
9: 0.90 cm**-1
10: 0.93 cm**-1
11: 0.95 cm**-1
12: 1.29 cm**-1
13: 1.30 cm**-1
14: 1.30 cm**-1
15: 1.37 cm**-1
...
29: 1.38 cm**-1
30: 1.40 cm**-1
31: 1.40 cm**-1
32: 1.40 cm**-1
33: 1.48 cm**-1
...
104: 1.48 cm**-1
105: 4.89 cm**-1
106: 4.95 cm**-1
107: 4.95 cm**-1
108: 5.03 cm**-1
109: 5.04 cm**-1
110: 5.04 cm**-1
111: 5.12 cm**-1
...
164: 5.12 cm**-1
165: 30.43 cm**-1
166: 36.35 cm**-1
167: 71.80 cm**-1
168: 227.53 cm**-1
169: 246.94 cm**-1
170: 265.41 cm**-1
171: 311.15 cm**-1
172: 418.07 cm**-1
173: 565.00 cm**-1
174: 608.49 cm**-1
175: 869.56 cm**-1
176: 901.82 cm**-1
177: 911.70 cm**-1
178: 1007.09 cm**-1
...

То есть первые 164 частоты в данном примере срежутся вообще любым критерием, и вся дискуссия сводится к учету трех частот 165-167

[Vit Nhoc]

Мне сказали, что интересующая меня информация есть в этой книге:

Глушко В. П. Термодинамические свойства индивидуальных веществ. 1 том.

Прошу подсказать, как её раздобыть или купить.

[Гесс]

Глушко В. П. Термодинамические свойства индивидуальных веществ. 1 том.
Прошу подсказать, как её раздобыть или купить.

В первом томе 3 книги. Проверьте имейл.

[alex7851]

Взгляните также на GoodVibes.

[Vit Nhoc]

Просьба ещё чуть подсказать. Я сделал расчёт энтропий, всё совпало с выдачей гауссиана. Также в энергию Гиббса вроде входит в каком-то виде теплоёмкость, и тут я запутался. Выше дали ссылку на сайт программы ASE, считающей термодинамические функции. По ней написано, что использовались формулы из книги C.J. Cramer. Essentials of Computational Chemistry, Second Edition. Wiley, 2004. В этой книге даны такие формулы:



Тут нет теплоёмкостей. На сайте же ASE приведены такие формулы:



Тут приведены даже не теплоёмкости, а их интегралы по температуре. В выдаче Гауссиана есть величины "CV":



По идее, эта CV должна быть интегралом из формулы ASE, только делённым на T. И у меня эти величины совпали для CV Translation и CV Rotational (формулы из сайта ASE и Гауссиан), но не совпали для CV Vibrational и вкладов в него по отдельным частотам. Что это означает? Может быть, Гауссиан выдаёт CV как простое произведение теплоёмкости на объем, и эта величина никак не входит в формулу для энергии Гиббса, поскольку в этой форме есть обозначенный интеграл? Если формула с интегралом как части суммы верная, как называется этот интеграл? Совсем забыл что учили по термодинамике, сорри.

[Гесс]

Это поможет?
https://www.youtube.com/watch?v=nh28_D4fIik
думаю для контроля можно использовать https://cccbdb.nist.gov/ihccalc1.asp - я очень надеюсь что в нисте закодировано правильно.
И да, поддерживаю прошлый ответ, https://pypi.org/project/goodvibes/ выглядит хорошо.

[YuraM]

Есть такая статья, thermochemistry in Gaussian, из серии Gaussian White papers или что-то такое, я все формулы оттуда кодил. Неплохо все дано в книге Molecular driving forces.

[Vit Nhoc]

Я наконец сделал этот расчёт термодинамических функций в Chemcraft, всё совпадает с выдачей других программ (по крайней мере с Gaussian). Провёл с помощью этого подхода расчёты pKa ряда карбоновых кислот и фенолов, и обнаружил, что для моего набора молекул оптимально установить threshold 300 cm-1 (если больше или меньше - согласие с экспериментом становится хуже). Интересный результат.
Просьба подсказать: на что сослаться в статье, когда я буду писать, что малые частоты дают большую ошибку в энергии Гиббса из-за ангармоничности, и где конкретно рекомендуется обрезать низкие частоты (или не обрезать а как это назвать по английски, может быть increase the frequency values to threshold value).

[Vit Nhoc]

Я сделал эти RRHO расчёты в Chemcraft, кому интересно предлагаю скачать последнюю версию:

https://chemcraftprog.com/download.html

На этой странице я выписал использованные формулы. Просьба взглянуть и посмотреть, вдруг вы там увидите ошибки:

https://chemcraftprog.com/help/thermodf ... lator.html

И ещё вопрос: в Орке колебательная энтропия считается не по стандартной формуле RRHO, а по формуле QRRHO (статья Grimme), можете подсказать в чём тут отличие?

[Гесс]

Вон есть еще mRRHO, как комбинация "complementary harmonic oscillator (S_v) and rigid-rotor (S_r) parts"
https://onlinelibrary.wiley.com/doi/10.1002/jcc.27129 (если нет доступа - пришлю)
Авторы ссылаются на уравнения 2-3 в https://pubs.rsc.org/en/content/article ... D1SC00621E

[kika]

Я хочу сделать расчёт термодинамических функций (энтропии, энергии Гиббса и пр.) молекул по молекулярным параметрам (структура, колебательные частоты), рассчитанным из квантовой химии.

Vit Nhoc, а Вы не могли бы мне посчитать энергию Гиббса для двух почти линейных молекул?
Или, может, Вы могли бы оценить, у какой из них эта энергия будет больше?
Они не очень сильно отличаются по структуре.
Заранее большое спасибо!

[Vit Nhoc]

Я хочу сделать расчёт термодинамических функций (энтропии, энергии Гиббса и пр.) молекул по молекулярным параметрам (структура, колебательные частоты), рассчитанным из квантовой химии.

Vit Nhoc, а Вы не могли бы мне посчитать энергию Гиббса для двух почти линейных молекул?
Или, может, Вы могли бы оценить, у какой из них эта энергия будет больше?
Они не очень сильно отличаются по структуре.
Заранее большое спасибо!

Вопрос немного странный - вам же выдаются эти энергии Гиббса в Гауссиане или Орке?
То, что я сделал в Chemcraft, нужно для другого - обрезать низкие или отрицательные частоты. Вам это не нужно.

[kika]

Вопрос немного странный - вам же выдаются эти энергии Гиббса в Гауссиане или Орке?
То, что я сделал в Chemcraft, нужно для другого - обрезать низкие или отрицательные частоты. Вам это не нужно.

Верно, если бы мне выдали эти энергии, то я бы не обращалась за помощью.
Видимо, это программные обеспечения, но я ими не располагаю.
Чтобы Вы посоветовали в этом случае?
Как грубо можно посчитать энергию Гиббса для органической молекулы?
Заранее спасибо!

[Vit Nhoc]

Вопрос немного странный - вам же выдаются эти энергии Гиббса в Гауссиане или Орке?
То, что я сделал в Chemcraft, нужно для другого - обрезать низкие или отрицательные частоты. Вам это не нужно.

Верно, если бы мне выдали эти энергии, то я бы не обращалась за помощью.
Видимо, это программные обеспечения, но я ими не располагаю.
Чтобы Вы посоветовали в этом случае?
Как грубо можно посчитать энергию Гиббса для органической молекулы?
Заранее спасибо!

1) Скачайте Гауссиан или Орку. На торрентах можно скачать Gaussian09, про Орку уже точно не помню, но она бесплатная;
2) Скачайте Chemcraft;
3) В каталоге Samples в Chemcraft найдите примеры входных файлов с типом задачи "оптимизация+частоты";
4) Нарисуйте в Chemcraft вашу молекулу (многие обычные молекулы есть в базе фрагментов через Edit/Add fragments), перейдите во вкладку Coords и скопируйте координаты атомов;
5) Скопированные координаты вставьте во входной файл;
6) Запустите этот расчёт. Можете использовать утилиту Tools/Scripts/Create bat file в Chemcraft, чтобы разобраться как в Гауссиане/Орке запускаются файлы;
7) Когда расчёт завершится, откройте выходной файл, там будет приведена энергия Гиббса молекулы.

[kika]

1) Скачайте Гауссиан или Орку. На торрентах можно скачать Gaussian09, про Орку уже точно не помню, но она бесплатная;
2) Скачайте Chemcraft;
3) В каталоге Samples в Chemcraft найдите примеры входных файлов с типом задачи "оптимизация+частоты";
4) Нарисуйте в Chemcraft вашу молекулу (многие обычные молекулы есть в базе фрагментов через Edit/Add fragments), перейдите во вкладку Coords и скопируйте координаты атомов;
5) Скопированные координаты вставьте во входной файл;
6) Запустите этот расчёт. Можете использовать утилиту Tools/Scripts/Create bat file в Chemcraft, чтобы разобраться как в Гауссиане/Орке запускаются файлы;
7) Когда расчёт завершится, откройте выходной файл, там будет приведена энергия Гиббса молекулы.

Vit Nhoc, спасибо большое!
Но я боюсь, что я это все не осилю - и установить, и разобраться. Но я рискну, если Вы будете рядом и поможете,
хотя, конечно, без Вашего терпения тут не обойтись.
1) вот тут можно скачать Gaussian09: https://gaussian-09w.software.informer.com/9.0/ ?
2) а Chemcraft тут скачивается: https://chemcraft.software.informer.com/1.8/ ?

[Vit Nhoc]

kika

Я вначале подумал что вы нашли пиратский сайт с ворованными копиями, но тут вроде надо вводить ключ. Как я понимаю, это какие-то мусорные сайты с файлопомойками, которые за счёт SEO отбирают у меня поисковый трафик и зарабатывают на рекламе. Гугл и Яндекс должны бороться с такими сайтами.
На этом сайте очень старая версия Chemcraft. Вы можете скачать последнюю версию моей программы по этой ссылке:

https://chemcraftprog.com/files/Chemcra ... _win64.zip

Ссылку для Gauusian09 не смотрел, можете попробовать, но подозреваю вам всё равно придётся идти на рутрекер. С ворованными программами есть такой аспект: вы можете посчитать что-то на ворованном Гауссиане, но публикация результатов без лицензии может иметь проблемы. Если же есть ворованный Chemcraft, то им можно пользоваться бесплатно, т.к. обычно его необязательно упоминать в статьях. Поэтому прошу соблюдать мои авторские права.