Альтернативная standalone программа для расчета сольватации

[Yury Minenkov]

Дорогие Друзья!

В группе квантово-химических расчетов ФИЦ ХФ РАН (https://vk.com/public221456203) была разработана альтернативная континуумная модель сольватации (solv). Модель позволяет вычислить энергию Гиббса сольватации произвольного растворенного вещества в произвольном растворителе, а также отдельные вклады в эту величину. Электростатический вклад в свободную энергию Гиббса сольватации вычисляется с использованием неитеративного приближения COSMO Vyboishchikov and Voityuk (DOI: 10.1002/jcc.26531) и фиксированных зарядов на атомах. Неэлектростатическая дисперсионно-отталкивательная энергия взаимодействия между растворенным веществом и растворителем рассчитывается через метод атом-атомных потенциалов Caillet−Claverie, использующий сеточный подход. Неэлектростатическая энергия Гиббса кавитации рассчитывается в рамках формализма Scaled Particle Theory (SPT). Радиус растворенного вещества может быть получен из трех схем: а) Pierotti−Claverie (PC), б) из площади молекулярной поверхности (SPT-S), в) из молекулярного объема (SPT-V). Радиус молекулы растворителя параметризуется через воспроизведение экспериментальных свободных энергий Гиббса сольватации 2530 нейтральных молекул в 92 растворителях. Предварительные тесты показывают, что подход SPT-V, основанный на зарядах, вычисленных в рамках модели CM5, дает наилучшие результаты. Метод может быть использован в качестве альтернативы моделям сольватации SMD (DOI: 10.1021/jp810292n) или uESE (DOI: 10.1002/jcc.26531).

Оригинальная публикация: https://doi.org/10.1021/acs.jctc.3c00410

Для запуска программы и получения результатов необходимы следующие шаги:

ШАГ 1: скачайте необходимые файлы:

1. Вспомогательные скрипты для создания input файлов для программы solv:

https://drive.google.com/file/d/1wSO1xf ... share_link

2. Сама программа solv:

https://drive.google.com/file/d/1zJvQox ... share_link

Разархивируйте папки (tar xvf ..)

ШАГ 2: Для запуска программы сольватации необходима x86_64 Linux система и входной input файл, который содержит Декартовы координаты Вашей системы и CM5 заряды на атомах в формате:

O -1.662892304 0.000000000 -3.660527232 -0.6615
H -1.662892304 0.759337000 -3.064484232 0.3307
H -1.662892304 -0.759337000 -3.064484232 0.3307

Файл в таком формате может быть получен с помощью любой квантово-химической программы, способной печатать Hirshfeld заряды. Мы ограничимся Priroda (Природа) и ORCA.

Priroda:

Создайте файл для вашей молекулы mol.inp, где включите в inp в группу $control:
$control
print=+charges
$end

Сделайте расчет в Priroda и получите файл mol.out
Далее, с помощью вспомогательных скриптов получите файл mol.HRS:

python2 ~/path_to_CM5/get_hrs_priroda.py mol.out

ORCA:
Создайте файл для вашей молекулы mol.inp, где включите в inp в группу %output:

%output
Print [P_Hirshfeld] 1
end

Сделайте расчет в ORCA и получите файл mol.out
Далее, с помощью вспомогательных скриптов получите файл mol.HRS:

python2 ~/path_to_CM5/get_hrs_orca.py mol.out

Полученный файл из ORCA или Priroda mol.HRS необходимо конвертировать в mol.M51, который будет содержать CM5 заряды.

python2 ~/path_to_CM5/hrs_to_m51.py mol.HRS

Скрипт использует оригинальную CM5 implementation (https://comp.chem.umn.edu/cm5pac/)

В конце у Вас должен появиться файл mol.M51 с CM5 зарядами.

ШАГ 3: Запустите программу сольватации на файле mol.M51:

~/path_to_SOLV/solv mol.M51 solvent_name grid1 grid2

solvent_name — название растворителя, см. solv.txt в папке SOLV для названия растворителей.

grid1 и grid2 — численные решетки для расчета. См. Оригинальную работу для деталей, согласно нашему опыту выбор grid1=300 and grid2 = 150 является оптимальным.

ПРИМЕР:
~/path_to_SOLV/solv mol.M51 toluene 300 150

РЕЗУЛЬТАТ:
Solvent: toluene // название растворителя
Rsolv = 2.873 (unscaled) // немасштабированный Rsolv из ван-дер-Ваальсова объема в Å
param: [0.882961, 0.964645, 0.997129] // коэффициенты для Rsolv в нашей параметризации для PC, SPT-S and SPT-V приближений
Solvent mass = 92.0626 // молекулярная масса растворителя в g/mole
Solvent den = 862.01 // плотность растворителя в kg/m3
Rsolv(vol) = 3.48574 // Rsolv из плотности растворителя, исключительно для сравнения
Gtr (1atm->x=1) = 3.21824 // поправка на термодинамические условия, при которых концентрация растворенного вещества в газе — 1 атмосфера, а в жидкости мольная доля x=1.
Warning: no R3 for: Sn, assign: X // предупреждения, если какой-либо параметр для какого-либо элемента отсутствует. В этом случае используется generic значение.
Warning: no k for: Sn, assign: X
// Все следующие значения для 1M стандартного состояния и в газе, и в жидкости.
Gel = -5.91 // Электростатическая часть энергии Гиббса сольватации в kcal/mol
Gdisp = -45.53 // Дисперсионная энергия в kcal/mol
Grep = 5.73 // Энергия отталкивания в kcal/mol
Rmol (SPT-S) = 6.24838 // радиус растворенного вещества из молекулярной поверхности в Å
Rmol (SPT-V) = 4.75975 // радиус растворенного вещества из молекулярного объема в Å
Gcav = 29.17 27.41 21.73 // энергия кавитации в приближениях PC, SPT-S and SPT-V в kcal/mol
Gsolv = -16.54 -18.3 -23.98 // Полная энергия Гиббса сольватации в kcal mol в приближениях PC, SPT-S and SPT-V. Мы рекомендуем значение SPT-V (-23.98 ккал/моль) как наиболее надежное.


ПРЕИМУЩЕСТВА нашей модели:
1) Легкая параметризация для практически любого растворителя. Необходимые параметры растворителя для параметризации: плотность, брутто-формула, давление насыщенных паров.
2) Исключена ошибка, обусловленная наличием зарядов вне полости (outlying/escaped charge)
3) Четкое разделение неэлектростатических вкладов на отдельные компоненты.

НЕДОСТАТКИ нашей модели:
1) На данный момент возможен только single-point расчет на фиксированной геометрии.
2) Не учитывается обратная поляризация растворенного вещества растворителем (Solute back polarization)
3) Специфические взаимодействия «растворитель — растворенное вещество» учитываются только неявно, через параметризацию радиуса растворенного вещества Rsolv.

Надеемся, модель будет полезна для Ваших исследований!

[amge]

Отлично! Работает.

Легкая параметризация для практически любого растворителя. Необходимые параметры растворителя для параметризации: плотность, брутто-формула, давление насыщенных паров.

А как это сделать? У меня много экспериментов в трифторметансульфокислоте (TfOH) и хотелось бы ее как растворитель для расчетов.
плотность 1.696 г/см³
брутто-формула CF3SO3H
давление насыщенных паров 1 мм рт. ст. при 25 гр.
Куда эти параметры вгонять?

На данный момент возможен только single-point расчет на фиксированной геометрии.

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

[Yury Minenkov]

Спасибо за вопрос.
Надо будет выложить код и объяснить процедуру параметризации. На сейчас просто добавьте в файл solv.txt:

name Triflic_Acid trifluoromethanesulfonic_acid TFMS TFSA HOTf TfOH
param 0.9851131439208984 1.0234050750732422 1.0449810028076172
den 1696 // AMGE (chemport)
elements S C H O F
index 1 1 1 3 3
eps 77.4 // https://arxiv.org/pdf/physics/9706028.pdf
diff 1e-9 // incorrect
vis 1e-9 // incorrect
rsolv 2.70156

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

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