Моделирование сольватации в воде большим количеством “экплицитных” молекул

[Vit Nhoc]

Я занимаюсь предсказанием pKa кислот в воде. Для этого надо рассчитать энергию нейтральной молекулы и энергию депротонированного аниона, разность между ними будет дескриптором pKa. Для более точных результатов нужно добавлять в модель “эксплицитные” молекулы воды (explicit water molecules). Обычно их добавляют 2-3, а я решил попробовать много, например 14 молекул:



Столкнулся с такой ситуацией: энергия очень зависит от того, какую конформацию найдёт расчёт (взаимного расположение молекул воды), и невозможно перебрать все конформации и найти наилучшую. Может быть правильно называть эти конформации геометрическими конфигурациями или как-то ещё? Мне кажется, хорошего результата теоретически можно добиться, если расчёты всех молекул проводить в одной и той же геометрической конфигурации, для этого необходима утилита в Chemcraft Structures merger.
14 молекул воды в данном случае – может быть это называется вторая сольватная сфера?
Пока мне не совсем понятно, насколько такой подход хорош; и надо очевидно искать другие статьи, в которых тоже авторы добавляли много молекул воды, причём для самых разных расчётов. Подскажите, как искать такие статьи.
Ещё вопрос: может быть для таких расчётов следует убрать имплицитную модель сольватации (PCM), т.к. с ней получится двойной учёт растворителя?

[bigM]

Я занимаюсь предсказанием...

извините, не удержался
https://www.youtube.com/watch?v=KfX3Tyf32bA

[Гесс]

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

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

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

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

Пока мне не совсем понятно, насколько такой подход хорош; и надо очевидно искать другие статьи, в которых тоже авторы добавляли много молекул воды, причём для самых разных расчётов. Подскажите, как искать такие статьи.

У Юрия Миненкова были статьи этого направления. Я был глубоко убежден что у Лайкова был скрипт/метод обкладывающий молекулу А заданным количеством молекеул Б с поиском оптимальной шубы, но я не нашел этой публикации. Amge посоветовал смотреть https://xtb-docs.readthedocs.io/en/latest/crestqcg.html, есть немало других утилит, например https://www.mpinat.mpg.de/628091/03-How_SOLVATE_works, https://chemrxiv.org/engage/api-gateway ... ethods.pdf

Ещё вопрос: может быть для таких расчётов следует убрать имплицитную модель сольватации (PCM), т.к. с ней получится двойной учёт растворителя?

нет, не получится.

[Vit Nhoc]

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

Сопутствующий вопрос: как всё-таки называть эти локальные минимум в англоязычной статье - local minimums, conformatins, configurations, geometrical configurations?

[YuraM]

PKa
Здравствуйте, опыта расчёта pka у меня нет. Из своего опыта, я бы предложил следующее. Запишите термодинамический цикл, в газовой фазе и растворе.
HAg = A-g + H+g, dGg
HAs = A-s + H+s, dGs
dGs=-RTlnKa
Ka = exp(-dGs/RT)
PKa = - lgKa
Т. Е. У Вас задача найти изменение dG реакции в растворе, dGs. Эта величина находится через сумму энергии Гиббса реакции в газе, dGg (только не забудьте термодинамическую поправку для протона и учесть разницу концентрации 1.89 kcal/mol) и разницу энергий Гиббса сольватаций продуктов и реактантов.
Энергию Гиббса сольватации протона берите из эксперимента, например, из справочника Маркуса. Энергию сольватации молекулы AH и иона A- считайте по кластерно-континуумной модели, лучше мономерный цикл. Кстати, для аниона не забудьте диффузных функции. Посмотрите как это делается в работах Pliego, или наши крайние работы. Мы предлагаем при использовании PCM ИЛИ SMD точное число молекул растворителя определять по термодинамический устойчивости комплексов молекул или ионов с растворителей.
A— + H2O = A—H2O, dG1
A—H2O + H2O = A—H2O2, dG2
Останавливайтесь, когда у Вас dGn будет больше 0. Мне кажется, Вы найдёте разное число молекул растворителя n для аниона A- и нейтральной молекулы HA. Рассчитайте энергии сольватации HA, A- и возьмёте эксперимент для протона, дальше найдёте PKa.
По-поводу накидать 14 молекул. Найти самую выгодную конфигурацию с 14 молекулами воды с HA и A- практически не реально, я бы даже dft dz энергия не верил, но даже они запретно дороги. А потом, почему 14? А не 100? Сколько молекул должно хватить, чтобы приблизиться к раствор? И, уверен на 100%, даже при 14 молекулах, самые стабильные конформации для AH и A-будут разными.

[YuraM]

Диффузные функции для всего нужны, но чувствую, что это будет ресурсоемкая задача...

[amge]

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

По моему, не стоит и пытаться. Даже если нашли расположение молекул воды, удачное для одного аниона, вовсе не факт, что такое расположение будет столь же удачным для другого аниона и тем более для нейтральной молекулы кислоты. Это будет вызывать большую неопределенность в разнице энергий. Например, вот кластер ацетат-аниона с 14 молекулами воды, полученный программой crest методом GFN2-xTB по алгоритму QCG.

AcO-14H2O.png

Сильно отличается от вашего. Попробуйте ради интереса заоптимизировать его и сравнить энергию с энергией своего кластера (координаты во вложении). Думаю, что будет большая разница, отнюдь не в доли ккал. А каждая ккал - это порядка единицы рКа.

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