Журнал неорганической химии, 2023, T. 68, № 11, стр. 1588-1598

Квантово-химические расчеты

Квантово-химические расчеты были выполнены в рамках метода теории функционала плотности по методике, принятой в наших предыдущих работах [29, 30]. Для вычислений использовали лицензионную компьютерную программу CRYSTAL17 [35, 36], в которой для построения базиса блоховских функций применяются гауссовы атомные орбитали. С целью передать взаимодействие между остовными и валентными электронами для атома S, а также остовными, валентными и субвалентными электронами для атома W использовали релятивистские псевдопотенциалы CRENBL [37, 38].

Для проведения расчетов был выбран гибридный обменно-корреляционный функционал HSE06 [39]. Сканирование зоны Бриллюэна проводили посредством сеток Монкхорста–Пака [40], состоящих из 18 × 18 × 10, 18 × 18 и 18 k-точек для объемного кристалла, монослоя и нанотрубок WS₂ соответственно. Одноэлектронные уравнения решали итерационно до самосогласования по энергии с точностью 3 × 10^–9 эВ.

Для воспроизведения ван-дер-ваальсовых взаимодействий между слоями в объемном кристалле, а также между стенками в многостенных нанотрубках в расчет включали дисперсионную поправку в приближении DFT-D2 [41].

Более детальное описание использованной методики приведено в работах [29, 30].

Выбор атомных потенциалов силового поля

Выбор функциональной зависимости для вкладов в потенциальную энергию является важнейшим этапом при создании силового поля. В данной работе мы остановились на потенциальной модели [31], предложенной нами ранее для моделирования многостенных нанотрубок на основе MoS₂. Эта модель содержит четыре типа потенциалов: 1) кулоновский, который выражается через эффективные заряды на атомах; 2) двухчастичный потенциал Морзе, отнесенный к короткодействующим взаимодействиям W–S и W–W; 3) потенциал Букингема, передающий парные взаимодействия между атомами серы; 4) трехчастичный потенциал Стиллинжера–Вебера, который обеспечивает вклад тройных взаимодействий, связанных с валентными углами WSW. Все упомянутые потенциалы реализованы в свободно распространяемой программе GULP [42], предназначенной как для молекулярно-механического и молекулярно-динамического моделирования, так и для калибровки параметров потенциалов. Конкретный вид использованных потенциалов указан в работе [31], a также представлен в электронном виде в дополнительных материалах. Найденные значения эмпирических параметров даны в табл. S1 . Там же результаты молекулярно-механических расчетов свойств калибровочных систем сопоставлены с их эталонными значениями (табл. S2 и S3 ). Данная потенциальная модель хорошо себя зарекомендовала [31] при моделировании структуры и свойств объемных фаз, монослоев, одностенных и многостенных нанотрубок на основе MoS₂.

Однако при тестировании указанной выше потенциальной модели выяснилось, что она неверно передает относительную энергию 2H- и 1T-фаз MoS₂, заметно завышая стабильность последней. Это обстоятельство приводит к неверным результатам молекулярно-динамического моделирования нанотрубок при повышенных температурах. Для устранения указанного недостатка в данную потенциальную модель был введен дополнительный трехчастичный потенциал, который понижает относительную выгодность октаэдрической координации атома металла по отношению к тригонально-пирамидальной, реально наблюдаемой в слоях WS₂. Данный потенциал заимствован из силового поля ESFF [43] и отражает угловую зависимость энергии для валентных углов, близких к 0°, 90° или 180°. В программе GULP [42] этот потенциал обозначен как Linear-three-body (LTB). В работе использовано следующее его представление:

где θ – валентный угол SWS, а k – эмпирическая постоянная.

Методика калибровки параметров

В настоящей работе для калибровки параметров силового поля использована библиотека EZFF [44]. Данная библиотека предназначена для многокритериальной оптимизации силовых полей, но, что самое важное, не привязана к конкретной молекулярно-механической программе, а может использовать GULP [42], RXMD [45] и т.д. Другим важным преимуществом этой библиотеки является то, что она распространяется свободно и имеет открытый исходный код, который, согласно лицензии, пользователь может редактировать. Данная библиотека написана на языке программирования Python3. Библиотека EZFF использует имплементацию генетических алгоритмов, реализованную в пакете Platypus [46].

Для визуализации результатов использовали библиотеки seaborn [47] и matplotlib [48]. Обработку и анализ полученных силовых полей проводили при помощи библиотек Pandas [49] и scikit-learn [50].

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

Вместо того, чтобы минимизировать погрешность каждой характеристики рассматриваемых структур отдельно, было решено разбить их на группы. Всего было выделено пять групп: 1) длины векторов трансляции решеток, 2) углы между векторами трансляции решеток, 3) характеристики упругости (модуль Юнга, элементы тензора упругости и т.д.), 4) частоты фононов, 5) относительные значения полной энергии. Для каждой из групп расчет погрешностей осуществлялся функцией MAPE (Mean Absolute Percentage Error) по формуле:

где $\widehat {{{y}_{i}}}$ – рассчитанное значение i-го параметра, ${{y}_{i}}$ – его истинное значение, полученное экспериментально или рассчитанное теоретически, а $\epsilon $ – произвольное небольшое, но строго положительное число, необходимое, чтобы избежать сингулярности, когда y_i = 0.

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

В настоящей работе процесс оптимизации силового поля разбивался на несколько этапов:

– сканирование пространства параметров с целью примерного определения областей, в которых параметры потенциалов оптимальны;

– анализ полученных результатов, выделение оптимальных силовых полей и разбиение их на группы по принципу наибольшей схожести;

– продолжение оптимизации параметров силовых полей внутри каждой группы с учетом большего числа модельных систем.

Из всех полей, полученных по ходу оптимизации, выбирали 1% полей с наилучшими значениями целевых функций. Выбранные поля подвергались финальному этапу оптимизации методом градиентного спуска в программе GULP [42].

Несмотря на то что алгоритм EZFF не может решить главную задачу – создание наборов готовых силовых полей, данный подход помогает найти хорошие начальные приближения, наличие которых заметно упрощает работу по созданию поля, пригодного для моделирования наносистем.

Калибровочные системы

Набор калибровочных (обучающих) систем определяется объектами предполагаемого моделирования – многостенными нанотрубками на основе WS₂. Очевидно, что большинство физических свойств нанотрубок потенциально содержится в свойствах исходных слоистых кристаллов, так как нанотрубки, по сути, воспроизводят ту же структуру, но в другой “свернутой” надмолекулярной геометрии. Поэтому в число эталонных структур были включены две стабильные фазы (2H и 3R) дисульфида вольфрама, а также его изолированные монослои. Кроме того, по указанным выше причинам к реальным фазам была добавлена гипотетическая 1T-фаза, в которой атом металла имеет октаэдрическую координацию. В набор эталонных свойств были включены постоянные решетки, координаты атомов, упругие постоянные, длины химических связей W–S и величины валентных углов SWS, WSW, а также (ограниченные) наборы фононных частот в различных точках зоны Бриллюэна. Для большинства указанных величин были использованы экспериментальные значения. В тех случаях, когда экспериментальные данные отсутствовали, использовали результаты квантово-химических расчетов (см. ниже).

Чтобы максимально корректно передать зависимость энергии взаимодействия стенок в многостенной трубке, в число калибровочных систем были добавлены шесть пар монослоев, расположенных друг от друга на разных расстояниях и с различным относительным сдвигом и ориентацией, что соответствует разной упаковке слоев в объемных фазах 2H и 3R [51]. Кроме того, были включены две нанотрубки с небольшими диаметрами: одна – типа “кресло” (6, 6), а другая – типа “зигзаг” (12, 0). Для указанных систем эталонными величинами выступали постоянные решетки и положения атомов, найденные квантово-химически. Для нанотрубок были дополнительно рассчитаны и включены первые шесть частот в центре зоны Бриллюэна для того, чтобы гарантировать отсутствие мнимых частот и правильно передать самые низкие частоты.

Для всех указанных выше калибровочных систем в число эталонных свойств входила их (потенциальная) энергия образования из объемной фазы 2H, рассчитанная квантово-химически.

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

Тестирование предложенного силового поля

Разработанное силовое поле было протестировано на одностенных и двустенных нанотрубках, которые предварительно были рассчитаны при помощи квантово-химического метода, описанного в предыдущем разделе.

Энергетическая стабильность нанотрубок традиционно оценивается по величине энергии сворачивания (E_str), которая вычисляется относительно энергии монослоя:

где E_NT и E_mono – полные энергии НТ и монослоя, рассчитанные на одну элементарную ячейку. Величины N_NT и N_mono (=1) составляют соответствующие числа формульных единиц.

Зависимости величин энергий сворачивания от средних диаметров нанотрубок, вычисленные посредством силового поля, сопоставлены с аналогичными зависимостями, полученными квантово-химическим методом, на рис. 2. При этом рассмотрены ахиральные нанотрубки типа “кресло” и “зигзаг” с диаметром в интервале 10–40 Å. Средний диаметр нанотрубок D_NT вычисляли как сумму (средних) радиусов цилиндров, соответствующих положениям центров атомов серы на внутренней и внешней оболочках. Как видно из рис. 2, согласие между молекулярно-механическими и квантово-химическими величинами достаточно хорошее: абсолютная ошибка составляет 3–5 кДж/моль.

Для подтверждения применимости силового поля при моделировании многостенных нанотрубок двустенные ахиральные нанотрубки с внешним диаметром в интервале 30–40 Å были рассчитаны как квантово-химически, так и в рамках силового поля. Оптимизировали геометрию двустенных трубок типа “кресло” (n₁, n₁)@(n₂, n₂) и “зигзаг” (n₁, 0)@(n₂, 0). Начальные структуры трубок создавали объединением двух одностенных нанотрубок одинакового типа хиральности и разных средних диаметров ${{D}_{{{{n}_{2}}}}}$ и ${{D}_{{{{n}_{1}}}}}$ с ограничением, чтобы исходный интервал между стенками ΔR_NT = (${{D}_{{{{n}_{2}}}}}$ – ${{D}_{{{{n}_{1}}}}}$)/2 составлял от 5 до 7 Å. Одновременно индексы хиральности n₁ и n₂ выбирали таким образом, чтобы при фиксированном межстенном расстоянии они имели максимально возможный общий делитель n. Такой выбор обеспечивает наибольшую симметрию двустенной трубки и, как следствие, снижает вычислительные затраты.

Как известно [29], одностенные ахиральные нанотрубки, свернутые из дисульфидов вольфрама, относятся к 4-му (типа “кресло”) или 8-му (типа “зигзаг”) семейству однопериодических групп, и у них отсутствует поворотная ось второго порядка, перпендикулярная оси трубки. Фактически это соответствует возможности двух различных ориентаций при упаковке монослоев в объемных фазах WS₂. Поэтому расчеты двустенных трубок проводили в двух вариантах: 3R-упаковка и 2H-упаковка, которые соответствуют либо одинаковой, либо противоположной ориентации трубок, свернутых из одного и того же монослоя.

Для сравнения стабильности трубок с различным межстенным расстоянием была рассчитана энергия связывания E_bind, которая показывает выигрыш в энергии (на одну формульную единицу) при объединении двух ОСНТ в одну ДСНТ. Эта величина отражает вклад межстенных взаимодействий в энергию двустенной трубки и определяется следующим образом:

где ${{E}_{{{\text{DWNT}}}}}$, ${{E}_{{{\text{SWN}}{{{\text{T}}}_{1}}}}}$, ${{E}_{{{\text{SWN}}{{{\text{T}}}_{2}}}}}$ – энергии ДСНТ и одностенных компонентов соответственно, отнесенные к одной трансляционной единице, а N_DWNT – число формульных единиц в элементарной ячейке ДСНТ. Помимо энергии связывания в число рассчитываемых величин для двустенных нанотрубок были включены энергия образования ДСНТ из объемного кристалла: которая, возможно, является самым практически значимым критерием стабильности нанотрубок, а также энергия образования ДСНТ относительно монослоя: которую можно непосредственно сравнить с энергией сворачивания E_str для одностенных нанотрубок близкого диаметра. В формуле (5) E_bulk – это рассчитанная энергия одной элементарной ячейки наиболее стабильной кристаллической 2H-фазы WS₂, содержащей N_bulk (= 2) формульных единиц.

В табл. 1 представлены перечисленные выше величины, полученные в результате квантово-химических и молекулярно-механических расчетов. Квантово-химические расчеты показали, что зависимость полной энергии двустенных нанотрубок от упаковки пренебрежимо мала (не превосходит 0.1 кДж/моль). Поэтому в табл. 1 представлены величины, полученные для трубок, оба компонента которых свернуты из монослоя в одной и той же ориентации. Согласие между квантово-химическими и молекулярно-механическими расчетами можно признать удовлетворительным. Для энергий образования различие между значениями, полученными с использованием обоих подходов, в среднем составляет 2 кДж/моль и не превышает 5 кДж/моль. Для энергий связывания средняя погрешность молекулярно-механических величин несколько больше (~3 кДж/моль).

Зависимость энергии связывания от межстенного расстояния показана на рис. 3. Как отмечалось ранее [28], исходное межстенное расстояние (расстояние между серединами стенок) ΔR_NT прямо пропорционально Δn_NT = n₂ – n₁. Разности индексов хиральности Δn_NT рассмотренных двустенных нанотрубок были подобраны таким образом, чтобы передать предполагаемый минимум на указанной зависимости. Сравнение эмпирических и квантово-химических результатов показывает их качественное согласие. Следует отметить, что при большом межстенном расстоянии оба подхода приводят к искаженной (сплюснутой) форме внешней нанотрубки с уменьшением межстенного расстояния на одних участках соприкасающихся поверхностей и с увеличением на других. По данным квантово-химических расчетов, оптимальное расстояние между стенками в двустенных трубках типа “кресло” составляет 6.4 Å, что соответствует Δn_NT = 7, 8, а это же расстояние в трубках типа “зигзаг” равно 6.2 Å, что примерно соответствует Δn_NT = 13. Молекулярно-механические расчеты дают следующие величины: 6.2 Å (Δn_NT = 7) для ДСНТ “кресло” и 6.0 Å (Δn_NT = 12, 13) для ДСНТ “зигзаг”.

Расчеты многостенных нанотрубок

В работе [25] представлено электронно-микроскопическое изображение высокого разрешения четырнадцатистенной нанотрубки WS₂ хиральности типа “зигзаг” с внешним диаметром 30 нм. Близкий диаметр имеют нанотрубки типа “зигзаг” с хиральностями около (300, 0). Для моделирования была выбрана серия нанотрубок с внешней стенкой хиральности (298, 0) и шагом хиральности Δn_NT = 12. Указанные параметры определяют следующую четырнадцатистенную нанотрубку (142, 0)@(154, 0)@(166, 0)@(178, 0)@(190, 0)@…@(286, 0)@(298, 0). На настоящем этапе работы выполнено моделирование внутренней части четырнадцатистенной НТ, образованной пятью одностенными нанотрубками от (142, 0) до (190, 0). Рассмотрена последовательность нанотрубок от двустенной (142, 0)@(154, 0) до пятистенной (142, 0)@(154, 0)@(166, 0)@(178, 0)@(190, 0). Наибольшая пятистенная трубка типа “зигзаг” содержит 4980 атомов в элементарной ячейке (рис. S1 ). Для сопоставления свойств нанотрубок с разным типом хиральности была рассчитана пятистенная НТ типа “кресло” с близким внешним диаметром и Δn_NT = 7: (82, 82)@(89, 89)@(96, 96)@(103, 103)@(110, 110). Элементарная ячейка этой нанотрубки содержит 2890 атомов (рис. S2 ). Как и в случае с пятистенной НТ типа “зигзаг”, был рассмотрен соответствующий набор одностенных и многостенных нанотрубок, связанный с пятистенной НТ типа “кресло”. Диаметры рассмотренных нанотрубок хиральности типа “кресло” близки по величине к диаметрам соответствующих нанотрубок хиральности типа “зигзаг”.

На рис. 4 показаны поперечные сечения трехстенных нанотрубок обоих типов хиральности с оптимизированными структурами. Видно, что НТ типа “зигзаг” имеет в сечении форму окружности, а НТ типа “кресло” приобретает граненую форму. Аналогичную форму имеют все стенки нанотрубок типа “кресло”, начиная с двустенных. Такое поведение многостенных нанотрубок наблюдалось ранее для целого ряда объектов, исследованных нами как квантово-химически, так и молекулярно-механически с использованием силовых полей. Граненая форма двустенных и трехстенных нанотрубок хиральности типа “кресло” на основе халькогенидов различных металлов (например, In₂STe [28] и SnS₂ [52]) обнаружена при проведении квантово-химических расчетов. Такой же результат показало моделирование с использованием силовых полей многостенных нанотрубок на основе MoS₂ [31], GaS, GaTe и Ga₂STe [34].

Зависимости энергий связывания (уравнение (4)), сворачивания (уравнение (6)) и образования (уравнение (5)) от числа стенок рассмотренных нанотрубок приведены на рис. 5. Приведенные графики показывают, что уже двустенные нанотрубки оказываются стабильными относительно монослоя. Настоящий результат согласуется с выводами, полученными ранее в работе [53], и подтверждает их. При этом многостенные нанотрубки WS₂ хиральности типа “зигзаг” стабильнее нанотрубок хиральности типа “кресло”, и относительная выгодность нанотрубок типа “зигзаг” нарастает с увеличением числа стенок. Энергия образования многостенных нанотрубок из объемной фазы положительна и с ростом числа стенок, естественно, стремится к нулю.

Важным структурным параметром многостенных нанотрубок является межстенное расстояние. Этот параметр экспериментально измеряется и может служить одним из критериев работоспособности силового поля и качества моделирования. В табл. 2 показаны межстенные расстояния для исследованных многостенных нанотрубок. По данным нашего моделирования, межстенные расстояния оказываются максимальными для самых тонких внутренних нанотрубок-стенок, уменьшаясь по мере роста их диаметров, но вновь возрастают во внешних частях многостенных нанотрубок. При этом все межстенные расстояния больше, чем расстояние между слоями объемного кристалла (6.162 Å) [51]. В работе [25] экспериментально обнаружено изменение межстенного расстояния в многостенной нанотрубке (30 стенок, внешний диаметр 500 Å) в зависимости от расположения стенок в середине, на внутреннем или внешнем крае нанотрубки. Однако приведенные значения [25] межстенных расстояний весьма велики: от 6.84 до 6.41 Å. По нашим данным, средние межстенные расстояния в нанотрубках типа “зигзаг” составляют 6.20 Å, в нанотрубках типа “кресло” – 6.25 Å. Наши величины почти совпадают с результатами измерений в работе [20], где в зависимости от метода синтеза нанотрубок их средние межстенные расстояния равны 6.214 Å (метод I) и 6.242 Å (метод II).

Диаметр одностенной НТ изменяется, когда НТ становится одной из стенок многостенной НТ. В табл. 3 приведены диаметры одностенных нанотрубок и их же диаметры в случае, когда эти НТ входят в состав многостенных нанотрубок. Согласно нашим расчетам, одностенные нанотрубки, входящие во внутреннюю часть многостенной нанотрубки, сжимаются. В середине многостенной нанотрубки изменения диаметров невелики, но по мере приближения к внешней поверхности диаметры стенок начинают превышать диаметры “свободных” НТ. Следует отметить, что эти изменения диаметров весьма малы и составляют десятые доли процента: 0.4–0.6%. Тем не менее современные экспериментальные методы позволяют измерять диаметр нанотрубок до десятых долей ангстрема (см. работы [22, 24]), и, следовательно, обсуждаемые изменения диаметров действительно могут быть обнаружены.