CC BY
102
В. А. Староверов
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОПТИМИЗАЦИИ ГЕОМЕТРИИ МОЛЕКУЛЯРНЫХ СИСТЕМ С ПОМОЩЬЮ ПРОГРАММНОГО КОМПЛЕКСА ACCELRYS
Ключевые слова: математическая модель, программный комплекс, компьютерное моделирование, молекулярная орбиталь.
В статье рассматривается математическая модель оптимизации геометрии исходных веществ и вычисления энергетического профиля реакции разложения с помощью программного комплекса accelrysmaterial studio 3.2, которая дает возможность быстрой прогностической оценки физико-химических свойств заданных веществ. Математическая модель дает возможность провести теоретические работы по компьютерному моделированию с использованием программно-аппаратного комплекса material studio 3.2. представляющая интерес в качестве разработке и использовании перспективных материалов при создании и поиске новых веществ, для изучения экологически чистых технологических добавок.
Keywords: mathematical model, software system, computer modeling, molecular orbital.
In this paper we consider a mathematical model of optimization of the geometry of the starting materials and calculate the energy profile of the decomposition reaction with the software package ACCELRYS MATERIAL STUDIO 3.2, which enables rapid prognostic assessment ofphysico-chemical properties of the specified compounds. The mathematical model makes it possible to carry out theoretical work on computer modeling using hardware-software complex MATERIAL STUDIO 3.2. of interest as the development and use of advanced materials in creating and finding new substances for the study of environmentally friendly processing aids.
Вычисления энергетического профиля реакции разложения возможно проводить с помощью программного комплекса accelrys material studio 3.2, который дает возможность быстрой прогностической оценки физико-химических свойств заданных веществ, представляющих интерес в качестве перспективных материалов для использования в создании экологически чистой технологии. Так, в плане поиска новых веществ, для изучения экологически чистых технологических добавок возможно провести теоретические работы по компьютерному моделированию с использованием программноаппаратного комплекса material studio 3.2:
1. Неэмпирические методы. Свободная энергия сольватации рассчитывается с учетом поверхности молекулы, доступной для атаки растворителя, то есть принимается во внимание конформа-ционное состояние молекулы.
2. Вычисление энергетического профиля реакции, поиск переходного состояния реакции и оптимизация геометрии переходного состояния в рамках полуэмпирических и неэмпирических методов:
- термодинамические вычисления в различных стандартных условиях (система с фиксированным числом частиц, объемом и энергией при изменяющейся температуре, число частиц, объем и температура системы постоянны, система обменивается энергией с окружающей средой, система с фиксированным числом частиц, давлением и энергией при изменяющемся объеме) в рамках полуэмпирических и неэмпирических методов.
- расчет энтальпии, энтропии и теплоемкости при различных значениях варьируемых переменных (температуре).
3. Расчеты функций Фукуи и вкладов атомов в граничные орбитали позволяют оценивать реакционную способность атомов в молекуле.
Для оптимизации геометрии исходных веществ в программе ACCELRYS MATERIAL STUDIO 3.2 использовался модуль DMol3.
Внутренняя схема оптимизации координат в модуле DMol3 обеспечивает хорошую оптимизацию для молекулярных систем. Данный модуль основывается на вычислении внутренних координат, способных скоординировать системы близко упакованных тел.
Внутренняя работа ратификации показывает, что это современное состояние внутренней схемы оптимизации координат периодических систем в 2-5 раз эффективнее, чем декартовские подходы, которые являются стандартом в вычислениях твердого тела на сегодняшний день.
В основу оптимизации молекулы лежит функциональная теория плотности, т.е. дискретное преобразование Фурье.
Функциональная теория плотности начинается с теоремы Hohenberg и Kohn, которая утверждает, что все свойства плотности описываются
полной энергией:
Е/М = тМ+Ф]+Е*сМ
(1.1)
где Т|р] - кинетическая энергия системы невзаимодействующих частиц плотности; и [р] - электростатическая энергия взаимодействия частиц плотности; Ехс[р] - обменно-корреляционная энергия.
Энергия корреляции, дискретное преобразование Фурье, явно зависит от всех количеств р.
При построении плотности нагрузки используются волновые функции, такие же, как и в других молекулярных орбитальных методах, но здесь волновая функция может использовать не симметричные молекулярные орбитали:
¥ = А(п)|ф(1)ф2(2)...фп (п)| (1.2)
Когда молекулярные орбитали ортонор-мальные:
(1.3)
Плотность возмущения определяется простой функцией:
2
где сумма Щ (r) i
РП=2
I
2
(1.4)
- это сумма по всем занятым
молекулярным орбиталям. Плотность, полученная от этого выражения также известна как плотность нагрузки. Молекулярная орбиталь может быть занята вращением наверх (альфа-электроны) или вращением вниз (бетта-электроны). Использование одного значения для альфа- и бетта-электронов вычисляется как ограниченное вращение. Использование же различных значений для альфа- и бетта-электронов вычисляется неограниченным вращением или поляризованным вращением. Во втором случае, можно сформировать две различные плотности нагрузки, одну для альфа молекулярных орбиталей, и одну для бетта молекулярных орбиталей. Их сумма дает полную плотность нагрузки, разность - плотность вращения, при избыточном количестве альфа и бетта вращений. Эти вычисления подразделяются на ограниченные и неограниченные [1].
Вычисление полной энергии
От волновой функции (1.1) и плотности возмущения (1.4) сумма энергий может быть записана как:
-V
2
nN - Z
и = 2 SUM----
(1.5)
(1.6)
+ - SS{фi (r )Ф j О2)-------------Фi 01)ф,- Os)) + 2 2
2 i J\ 1 J 2 r - r 1 J ! aR<,
ap<a\Ra -Rf}\
N
-a(p(i)L
- {- р(г1)гы} +(р(г1 + гж
В формуле (1.6) 2а - относится к нагрузке на ядре системы N - атомов. Первый член гУМ относится к электронному ядру. Второй гУМ представляет вращение электронного облака. Заключительный элемент УММ описывает вращение ядра.
Вычисление аппроксимации обменно-
корреляционной энергии
Расчет полной энергии (1.1) базируется на обменно-корреляционной энергии, которая требует некоторой аппроксимации для этого метода и должна быть достаточно универсальной, чтобы обрабатывать на вычислительной технике. Достаточно простой и хороший метод аппроксимации - локальное аппроксимирование плотности, которое базируется на известном методе обменно-корреляционной энергии однородного электронного газа. Аналитические исследования были сделаны несколькими
исследователями (Hedin, Lundgvist, Ceprrley, Adler, Barth). Локальная аппроксимация плотности предполагает, что плотность нагрузки изменяется медленно в масштабе атома. Полная энергия обменной корреляции может быть получена объединением однородных газовых облаков и описывается уравнением:
sxc [p ]= J p(r )e xc
\_p{r )] (1.7)
где Sxc [p ] - энергия обменной корреляции частицы в однородном электронном газе, r - число частиц [2;3].
Общая плотность вращения
Самая простая форма потенциала обменной корреляции - форма, полученная Slater (1951), который использовал Sxc [p] =р1/3. в этом приближении, корреляция не включена. Более сложные приближения используют формы, полученные Vosko, и другие, которые затруднительно использовать в вычислительной технике [1].
Расширение градиента плотности
Следующим шагом в улучшении локальной плотности вращения (LCD) является модель, учитывающая неоднородность электронного газа, который естественно находится в любой молекулярной системе. Это может быть достигнуто расширением градиента плотности, и относится как нелокальное приближение плотности вращения (NLSD). За последнее время данный метод хорошо себя зарекомендовал и подтвердил, что градиент корректировал обменно - корреляционную энергию, Exc[p,d(p)/dr], необходимую для изучения термохимии молекулярных процессов.
Нелокальное приближение плотности вращения (NLSD) функции используют при обобщенном приближении градиента для корреляционной функции.
Вычисление полного выражения энергии
С учетом вышесказанного полная энергия может быть записана как:
(1.8)
Et Ы=5(ф,- +
ф ) + (p(ri) Sxc[p<>1) + —р-VN] |) + VNN
Чтобы определить фактическую энергию, изменения в Е1 должны быть оптимизированы относительно изменений г, и приведены к ортонормаль-ным ограничениям [4]:
SEt
'-22 s..
Sp i j
!JV /
= 0
(1.9)
где Sy - лагранжевые множители.
Уравнения Kohn-Sham
Этот процесс ведет к набору соединенных уравнений, первыми предложили Kohn и Sham.
|-V
2
2
-Vn + Ve +<%c[ p]m = siфi
(1.10)
где дхс - потенциал обменной корреляции, который следует из дифференциала Ехс. Для приближения
2
R
R
2
локальную плотность вращения аппроксимируют, где потенциал ^ можно записать как:
Мхе = “ (р£хе) (1.11)
Использование собственных значений для быстрого нахождения энергии из уравнения (1.10) можно представить как [5]:
(1.12)
Е = Х ег + уР
Вычисление расширения молекулярных орбиталей в терминах основных функций
На практике очень удобно вычислять молекулярные орбитали с помощью атомов.
К-=2С,
(1.13)
Атомные орбитали Хд называют основными атомными функциями, где Сщ - коэффициенты расширения молекулярных орбиталей. Возможно несколько вариантов для основного набора, включая Гауссовские функции, волновые и числовые орбитали используемые при математическом моделировании. В последних расчетах при математическом моделировании используются волновые и числовые орбитали.
В отличие от молекулярных и атомных орбиталей - не ортонормальных, уравнение (1.10) можно записать в виде:
ИС=е8С, (1.14)
где
(1.15)
/ 1 2 \
№(r1) —2 vn + +Ve + MxcPr\ Xv (r1) )
S цу = {x^(r\)\ Xv(r\)) (1.16)
Оценка выражения (1.15) должна быть достигнута в соответствии с 3D числовой интеграцией и с учетом потенциала обменной корреляции. Данные элементы должны быть приближены к конечным суммам:
нцу=2 Хц(п )HeJf(r )У )w(ri) (1.17)
Это является дополнительной информацией для оценки данного выражения в соответствующих модулях программы Accelrys.
Способность оценивать производную полной энергии относительно геометрических изменений является актуальным для изучения химических элементов. Без данных исследований процедура вычислений становится тяжелой и очень трудоемкой, требующая привлечения большого количества квалифицированных специалистов. Данная теория является успешным методом для прогнозирования химических структур [6].
Эффект сольватации COSMO
Эффект COSMO позволяет исследовать эффекты сольватации.
На рисунке 1.1 представлена модель сольватации, в которой молекула растворенного вещества образует полость внутри диэлектрического кон-
тинуума диэлектрической постоянной, е, представляющего растворитель. В результате распределения зарядов растворенного вещества диэлектрическая среда поляризуется. В отличие от других анализов с помощью CSM COSMO не требует решения довольно сложных граничных условий диэлектрика для получения зарядов поляризации. Эти заряды рассчитываются с использованием значительно упрощенного граничного условия проводника. Заряды пересчитываются с учетом коэффициента f(e)=(e-1)/(е + 1/2) для получения приемлемой аппроксимации для зарядов поляризации в диэлектрической среде.
Поляризуемый растворитель
Заряды поляризации
Рис. 1.1 - Модель сольватации
Отклонение при точном решении приближенного COSMO являются незначительными. Для сильных диэлектриков вроде воды это отклонение меньше чем 1 %, в то время как для неполярных растворителей с е~2 оно может достигать 10% от общего эффекта поляризации. Однако для слабых диэлектриков, этот эффект незначителен, поэтому абсолютная ошибка составляет меньше чем 1 ккал/моль.
В целом COSMO - значительно упрощает подход центрального сервисного модуля (CSM) без снижения точности. В связи с таким упрощением COSMO позволяет более эффективно применять CSM в квантово- химических программах и точнее вычислять градиенты, которые позволяют оптимизировать геометрию растворенного вещества в пределах диэлектрической проницаемости.
Заряды поляризации определяются по граничному условию потенциала исчезновения на поверхности проводника. Если принять q за вектор зарядов поляризации на поверхности полости, а Q=p+Z - за общие заряды растворенного вещества, например, плотность электронов, р, и ядерные заряды, Z, то вектор потенциалов Vtot, на поверхности выражается в виде:
Vtot=BQ+Aq=Vsol+Vpol, где BQ - потенциал, обуславливаемый зарядами растворенного вещества; Q, и Aq - потенциал, обуславливаемый зарядами на поверхности, q; B и A - матрицы Кулона. В случае проводника приемлемо уравнение Vtot=0, при этом заряды поляризации определяются как:
q=-A-1BQ (1.18)
COSMO позволяет исследовать свободную энергию сольватации методами электростатики.
и
Кроме того, на общую свободную энергию сольватации влияют неэлектростатические процессы; они описывают взаимодействие дисперсий и эффект образования полости [7].
Термодинамические вычисления
Результаты термодинамического анализа могут использоваться для вычисления энтальпии образования (И), энтропии (8), энергии Гиббса (в) и способны проводить вычисления при высоких температурах в постоянном давлении в изотермическом процессе (ср), где полная энергия приводит к полной электронной энергии при 0 К. Для вычисления И,8,в и ср при начальных и критических температурах используются различные методы вычисления, которые будут описаны ниже.
Формулы базируются на работе И1гапо. В каждом случае описаны два выражения для переменных вкладов: одно для линейных и одно для нелинейных систем.
Расчет энтальпии
Расчет энтальпии в идеальном газе рассчитывают так:
И(Т)=БУ1Ь(Т)+ЕГо4(Т)+Б4гаП8(Т)+КТ, (1.19)
где Я - универсальная газовая постоянная.
В стандартных условиях энтальпия описывается соответственно для вибрационных, вращательных и переводных вкладов. Вклады имеют следующие значения:
3
C
tra
R
Etra = 2 RT
Erot (linear ) = RT
3
(1.20)
Erot(nonlinear ) = — RT
2
Evi
R 1 h + R hVi
k 2 “ ' k , [l - exp(-hVj / kT)]
R^ hvj exp(-hvt / kT)
где к - постоянная Больцмана, И - постоянная Планка, у1 - индивидуальные вибрационные частоты.
Расчет энтропии
Вклад в энтропию рассчитывают по форму-
лам:
Stra = 5RlnT + 2Rlnw - Rlnp - 2.3142S Srot (linear) = R lnl S'~ 2'T I + К
Srot (nonlinear) = — ln
(1.21)
+ і К 2
2 -/а И И И ^ Не)
Иу, / кТ ехр(-Иу, / кТ) г т
и = ЯI—--------1--------- - ЯЕ 1п|1-ехр(-Иу, / кТ)1
у,Ь , 1 - ехр(-Иу / кТ) , 1 ^ г п
где w - молекулярный вес, 1х - момент инерции на оси х, с - число симметрии.
Расчет теплоемкости.
Вклад теплоемкости при постоянном давлении описывается формулами, приведенными ниже:
2
(1.22)
Crot (linear) = К
3
Сг-ot (nonlinear) = — К
С = (hvi / kT)2exp(~hvi /kT)
Vlh i [l - exp(-hvi /kT)]
Для выгчисления АE реакции необходимо сложить полную энергию реагентов и энергии компонентов. Для выгчисления АH(T) необходимо добавить к выгаисленному теплосодержанию электронную энергию каждого компонента [8-9].
При значениях энтальпии и энтропии, рассчитанных по универсальной химической программе accelrys material studio 3.2, при заданных температурах легко рассчитать энергию Гиббса.
А G = А H - T А S (1.23)
По формуле 1.24 можно рассчитать энергию реакции при различной диэлектрической проницаемости:
реакции- продуктов 'исход, (1.24)
где АGреаKцИИ - энергия Гиббса реакции, ккал/моль, EАGпPOдVKTOв - энергия Гиббса продуктов реакции, ккал/моль, ^G^^ - энергия Гиббса исходных веществ, ккал/моль.АGреакции характеризует возможность самопроизвольного протекания реакции.
По формуле 1.25 можно рассчитать энергию переноса при различной диэлектрической проницаемости.
^перенося средьг ^GBgKyyM{ (125)
где АОпереноса - энергия переноса Гиббса при различной диэлектрической проницаемости, ккал/моль, АGсредЫ - энергия Гиббса в среде, ккал/моль, АGвакуум
- энергия Гиббса в вакууме, ккал/моль, АGперенOCa характеризует легкость протекания процесса.
Выводы
1. Установлено, что такие характеристики как АS, АИ, АG, АGреакции и АGпереноса позволяют достоверно оценить глубину протекающей реакции, энергозатраты, длительность и глубину протекающего процесса взаимодействия компонентов при обработке порохо-вык элементов.
2. С помощью универсальной химической моделирующей программы ACCELRYS material studio 3.2. можно провести расчеты по определению термодинамических характеристик исследуемык веществ.
3. С помощью универсальной химической моделирующей программы ACCELRYS material studio 3.2. можно оптимизировать геометрию молекул исследуемых веществ.
4. При нахождении энергии Гиббса можно рассчитать энергию реакции при различной диэлектрической проницаемости и рассчитать энергию переноса при различной диэлектрической проницаемости.
П Sk cIa 8л cln Sn-2clc (kT
Литература
1. Boese, A.D., Handy, N.C.J. Chem. Phys., 114, 254 (2001)
2. Ziegler, T. “Approximate density functional theory as a practical tool in molecular energetics and dynamics”, Chem.Rev., 91,451 (1991)
3. Hedin, L. Phys. Rev. A, 139, 350 (1965)
4. Delley, B. “Analytic energy derivatives in the numerical local-density-functional approach”, J. Chem. Phys., 94, 158 (1991)
5. Kohn, W.; Sham, L.J. “Self-consistent equations including exchange and correlation effects”, Phys. Rev. A, 140,260 (1965)
6. Bergner, A.; Dolg, M.; Kuechle, W.; Stoll, H.; Preuss, H. Mol. Phys., 80, 351 (1993)
7. Bayly, C.I.; Cieplak, P.; Cornell, W.D.; Kollman, P.A. “A well-behaved electrostatic potentional based method using charge restraints for deriving atomic charges: the RESP model”, J. Phys. Chem., 97, 215 (1993)
8. Pulay, P. Chem. Phys. Lett., 73, 393 (1980)
9. Косточко А.В., Проблемные вопросы утилизации поро-хов и некоторые области применения их в народном хозяйстве Косточко А.В., Косточко А.А., Ибрагимов Р.А., Храмова Е.В. «Вестник Казанского технологического университета» Казань 2012 том 15 №10 с 54
© В. А. Староверов - д.т.н, профессор кафедры ХТВМС КНИТУ, biserton@mail.ru.
