CC BY
27
УДК 517.958:518.12(075.8)
А.М. Попов, Г.Н. Шумкин2
МНОГОМАСШТАБНАЯ МОДЕЛЬ МОЛЕКУЛЯРНОГО ПЕРЕКЛЮЧАТЕЛЯ*
В работе представлено математическое моделирование "из первых принципов" молекулярного переключения на основе индуцированной током сканирующего туннельного микроскопа реакции изомеризации в молекуле нафталоцианина. Поверхность свободной энергии Гиббса и путь реакции анализируются методом метадинамики в рамках численного квантово-механического кода молекулярной динамики Кар-Парринелло CPMD. Расчеты проведены на супер-ЭВМ IBM Blue Gene/P, установленной на факультете ВМК МГУ. Вдоль введенной координационной переменной найдена высота энергетического барьера, который следует преодолеть для прохождения реакции изомеризации. Предложена многомасштабная модель молекулярной динамики Эрнфеста на основе возбужденных электронных состояний. Модель использует координационное направление реакции на поверхности свободной энергии, полученной в расчетах полной молекулы.
Ключевые слова: молекулярные переключатели, квантовая молекулярная динамика, вычислительные многомасштабные модели, путь химической реакции.
1. Введение. Вычисления "из первых принципов" широко используются для моделирования молекулярных переключателей [1-3]. Создание молекулярных электронных устройств, которые можно переключать между состояниями с высоким и низким сопротивлением, является одним из важнейших направлений молекулярной электроники. Так как миниатюризация этих технологий продолжает развиваться, то фундаментальная проблема идентифицирования и понимания наименьших физических систем, которые способны к переключательному поведению, вызывает растущий интерес [1]. Одиночные молекулы, несущие электрический ток между парой металлических электродов, могут быть использованы как электронные компоненты, такие, как логические затворы, элементы памяти, транзисторы или переключатели [1, 2].
Настоящая работа посвящена математическому моделированию "из первых принципов" экспериментов [3] по молекулярному переключению проводимости в молекуле нафталоцианина под действием тока сканирующего туннельного микроскопа (STM). Эксперименты показывают, что процесс молекулярной изомеризации может быть использован в рамках молекулы без изменения ее молекулярной геометрии. Бистабильность в позиции двух водородных атомов во внутренней области одиночной молекулы нафталоцианина задает двухуровневую систему.
Для компьютерного моделирования "из первых принципов" процесса изомеризации в свободной молекуле нафталоцианина мы используем молекулярную динамику Кар-Парринелло [4] и точный гамильтониан молекулы. Этот метод реализован в численном квантово-механическом коде молекулярной динамики CPMD (Car-Parrinello Molecular Dynamics) [4, 5]. Движение всех ядер используется для нахождения пути реакции. Основное внимание уделяется вычислению барьера на поверхности свободной
1 Факультет ВМК МГУ, проф., д.ф.-м.н., e-mail: popovQcs.msu.su
2 Факультет ВМК МГУ, асп., e-mail: georgiy-shQyandex.ru
* Работа выполнена при поддержке гранта РФФИ 09-07-12068.
6 ВМУ, вычислительная математика и кибернетика, № 3
энергии (FES), который необходимо преодолеть для прохождения реакции переключения. Нахождение пути реакции на рифленом ландшафте свободной энергии представляет важный шаг в понимании механизма реакции изомеризации. Поиск пути химической реакции представляет собой трудную задачу нахождения седловых точек в многомерном пространстве.
Для такого поиска мы используем метод метадинамики. Это мощный метод для эффективного вычисления поверхности свободной энергии как функции ограниченного набора коллективных переменных с использованием специальной стохастической молекулярной динамики [6]. Метод метадинамики реализован в CPMD и используется нами для анализа поверхности свободной энергии Гиббса.
Вычисления проводились на суперкомпьютере IBM Blue Gene/P (установленном на факультете ВМК МГУ), и продолжительность каждого варианта расчета составляла порядка 20 ч. Для оптимизации молекулярного переключателя требуется проведение многих расчетов. Поэтому необходимость предложенной в работе упрощенной, но многомасштабной модели очевидна.
2. Основные уравнения и постановка задачи. При построении моделей молекулярной динамики "из первых принципов" исходят из уравнения Шредингера и кулоновского взаимодействия заряженных частиц. Мы используем математическую модель молекулярной динамики Кар-Парринелло [4] на основе приближения Борна-Оппенгеймера. В приближении Борна-Оппенгеймера считается, что движение массивных ядер может быть описано в классическом приближении. Силы, действующие на ядра со стороны электронов, рассчитываются на основе решения уравнения Шредингера для электронов в основном состоянии. Мы будем использовать квантово-механическую формулировку, основанную на уравнениях функционала плотности (БРТ [7]). Уравнения модели Кар-Парринелло записываются в виде
MjRJÍÍ) =
5
д
■ЕКЗ[Ш,{Щ}}
5фг
ÖRj
EKS[{ipi},{ +
з i
(1)
где М/ и К/ — масса и координата ядра I, р, — фиктивная масса электрона, — орбитали Кона-Шэма, /-.'д с — энергия Кона-Шэма. Последний член в уравнении (1) характеризует голономные связи, которые имеют место в связи с требованием ортогональности орбиталей. Это приводит к необходимости на каждом шаге моделирования вычислять действующие на ионы силы (— -¡^-¡Екв)-, действующие на орбитали силы (—-¡^Екв) и множители Лагранжа Лц.
Энергия Кона-Шэма определяется выражением
ЕКзшл ад = Ёл / Гг(г)
где Uext — ионный потенциал и
dr+
Uext(r)p(r)dr+l f f drdr1 + Exc[p] + Ej[{Rf}],
r — r
— плотность заряда электронов, — число заселенности, Ехс — обменно-корреляционная энергия, /•.'/ — электростатическое ион-ионное взаимодействие.
Используется разложение волновых функций по плоским волнам. Орбитали Кона-Шэма записываются в виде
1
.¡Gr
G
где — объем ячейки, в — вектор обратной решетки, сДв) — фурье-коэффициенты орбитали г. В СРМБ используется быстрое преобразование Фурье для преобразования величин между двумя пространствами.
3. Моделирование поверхности свободной энергии "из первых принципов" методом метадинамики. Метод позволяет нам найти профиль FES вдоль координационной переменной и получить энергетический барьер реакции изомеризации. В то же время это позволяет изучить эволюцию электронной плотности в процессе метадинамики. Оптимизированная геометрия молекулы нафтало-цианина в состоянии первого изомера показана рис. 1, a, a в состоянии второго изомера — на рис. 1, б. После стадии термолизации мы находим термодинамическое равновесие, которое соответствует определенной температуре.
Рис. 1. Оптимизированная геометрия молекулы нафталоцианина на различных стадиях процесса метадинамики, вычисленная с помощью кода СРМБ: начальная конфигурация молекулы нафталоцианина (а); конечная конфигурация, обладающая геометрией второго изомера (б); рассчитанная структура молекулярных орбиталей Ы1МО (в) и Ы1МО + 1 (г). Белые шары соответствуют ядрам водорода, черные — ядрам азота, серые — ядрам углерода
Путь химической реакции обычно анализируется в терминах координационных переменных. Внутренняя координата реакции (IRC) [8] определяется как наикратчайший путь на поверхности потенциальной энергии, определяемой массами ионов, который соединяет состояние химического реагента с состоянием продукта реакции. Для изучения реакции водородной изомеризации, происходящей в молекуле нафталоцианина, была предложена координационная переменная, которая принимает во внимание координаты двух атомов водорода и двух атомов азота в центральной полости молекулы (рис. 1, а). Определим вспомогательную переменную в следующей форме:
8l(HuH2,NuN2) =
1 -
^RH1N1 ^ 6 ^ _ ^RH1N2 J _ ^RH2N1 J _ ^RH2N2 ^
_ ( Rh1N1 \ 12 _ ( Rh±N2 \ 12 -j _ f Rh2N± \ 12 -j _ f Rh2n2 \ 12
V Ro J 1 { Ro J 1 { Ro J 1 { Ro J
> .
Здесь RhiNj — расстояния между г-м атомом водорода и j-м атомом азота. Расстояние, задаваемое Rq, есть равновесное расстояние между атомом водорода и атомом азота, равное 1.05 Ä. Такая координационная переменная является несимметричной по углу вращения атома водорода. Поэтому была
использована координационная переменная в следующей форме:
s2 = s1(H1,H2,N1,N2) - s1(H1,H2,N3,N4).
В процессе молекулярной динамики вычисляется IRC для получения переходного состояния. Затем эволюция реакционной системы должна быть исследована вдоль минимального пути реакционной системы. Целью метода метадинамики является эффективный анализ поверхности свободной энергии F(s) как функции ограниченного набора коллективных переменных sa (длин связей и торсионных углов). Свободная энергия записывается в следующей форме:
где кв — константа Больцмана, Т — абсолютная температура и Р — функция распределения:
где и (г) — потенциал взаимодействия. Вводится фиктивная частица с координатой ,за для каждой коллективной переменной и используется техника расширенного лагранжиана. Теория функционала плотности (БРТ) используется для вычисления электронной структуры с использованием молекулярной динамики Кар-Парринелло [4] (СРМБ). В полной кинетической энергии фиктивная частица с достаточно большой массой р,а адиабатически отделена от ионных и электронных степеней свободы. Фиктивная частица ,за связана с ее действительной коллективной переменной 5(г) гармонической силой. Искусственный зависящий от времени потенциал ^(в, I) добавляется в систему вдоль траектории:
где вг = {«^(¿г)}- Этот потенциал есть сумма отталкивающихся потенциальных холмов в форме гаус-сианов, каждый из которых имеет высоту Н. Метадинамика ,за определяет эффективность, с которой исследуется поверхность свободной энергии. Затем потенциальные холмы добавляются в зависящий от истории потенциал ^(в, I) с меташагом тмтв, который на один или два порядка больше, чем шаг молекулярной динамики тмв- Потенциал вынуждает систему посещать те точки в конфигурационном пространстве, которые соответствуют яме реагента. Метадинамика искусственной частицы добавляет искусственный потенциал (холмы) в состояние реагента до тех пор, пока этот потенциал не скомпенсирует первую яму свободной энергии. Теперь система может найти наинизшее переходное состояние к следующему локальному минимуму (продукту реакции). Когда все локальные минимумы свободной энергии заполнены холмами, система может свободно двигаться от состояния реагента к состоянию продукта. В результате профиль -Р(в) может быть получен с произвольной точностью как взятый с обратным знаком потенциал ^(в, I).
Вычисленный профиль поверхности свободной энергии вдоль координационной переменной показан на рис. 2. Сходимость процедуры достигнута специальным набором параметров метадинамики. В действительности эта сходимость была получена проведением последовательности различных расчетов, каждый из которых стартовал с ситуации, имевшей место на предыдущем расчете с предыдущими параметрами.
Два равновесных состояния соответствуют величинам координационной переменной « = 0.6 и 5 = —0.6. Значение 5 = 0.0 относится к переходному состоянию. Вычисления состояли из четырех последовательных вариантов с различными параметрами. Первый расчет был проведен с высотой гауссиана Н = ЗквТ и полушириной Ш = 0.03. Заметим, что в этих вычислениях полуширина не изменялась и температура была равна Т = 30 К. Конец первого расчета соответствовал состоянию, когда система достигала второго равновесного состояния. Затем поверх вычисленного потенциала Ус добавлялись гауссианы с в два раза меньшей высотой Н = 1.5квТ, затем опять проводилась термоли-зация системы. Таким образом, набор вычислений был выполнен с высотами гауссианов Н = 0.75квТ и Н = 0.375квТ. Шаг метадинамики был равен 800 шагам молекулярной динамики, шаг тмв был 0.1125 Гб. Полученная высота барьера равна 52 Ы/то1 с точностью вычислений, равной 1 Ы/то1.
Видно, что зависимость свободной энергии имеет два основных минимума и соответствует двум состояниям переключения. Барьер между этими двумя состояниями должен быть преодолен за счет зарядки молекулы или добавления зарядов на Ы1МО или как действие электрического поля. Влияние силы в направлении координационной переменной приводит к вытягиванию водородной связи
F(s) = ^kBT\n[P(s)}
EES, kJ/mol
-60 -1-'-'-1-'-'-1-
-0,8 -0,6 -0,4 -0,2 0 0,2 0,4 0,6 5
Рис. 2. Вычисленный профиль поверхности свободной энергии вдоль координационной переменной 5. Сходимость достигнута в результате четырех последовательных расчетов с параметрами: 1) Н = 3квТ, Дз = 0.03; 2) Н = 1.5/свТ; 3) Н = 0.75квТ; 4) Н = 0.375/свТ
и торсионных углов в начальном равновесии. Когда наступает переход через барьер, старые связи обрываются и образуются новые связи. Этот процесс соответствует переключению молекулы.
Внутреннее вращение в молекуле есть важная часть процесса изомеризации. Эффект конечной температуры сказывается на колебаниях ядер в молекуле и изменяет профиль поверхности свободной энергии.
4. Редуцированная модель переключения на основе возбужденных молекулярных состояний. Высота энергетического барьера и возможность его преодоления (с помощью электрических и магнитных полей) определяется составом молекулы и связями атомов в ней. Возможность вычисления высоты барьера определяет оценку изменений в композиции молекулы для удовлетворения технологическим требованиям. Вычисления проводились на суперкомпьютере IBM Blue Gene/P (установленном на факультете ВМК МГУ) для 82-280 атомов в молекуле, и длительность одного варианта вычислений для определения барьера реакции составляла 20 ч. Для оптимизации молекулярного переключателя требуется проведение многих расчетов. Поэтому необходимость разработки упрощенной, но многомасштабной модели очевидна. Мы исходим из того, что основная трудоемкость вычислений требуется для определения поверхности свободной энергии и барьера реакции для полной молекулы. Если эта часть работы выполнена с помощью CPMD, то быструю часть переключательного процесса можно построить на основе упрощенной динамической модели молекулярной динамики Эрнфеста. Упрощение состоит в том, что основная конфигурация рассчитана CPMD и не меняется в быстром процессе молекулярной динамики Эрнфеста [9]. Более того, предлагается использовать только несколько наиболее важных, возбужденных орбиталей Кона-Шэма, показанных на рис. 1, в, г, где представлена рассчитанная структура молекулярных орбиталей LUMO (рис. 1, в) и LUMO + 1 (рис. 1,г). Видно, что структуры незанятых орбиталей LUMO и LUMO + 1 повернуты на 90° по отношению друг к другу.
Видно, что процесс переключения может быть связан с возбужденными орбиталями следующим образом. Действие внешнего потенциала за счет STM приводит к возбуждению тока через молекулу и электронного потока через LUMO- и LUMO + 1-орбитали, создавая возмущение электронной плотности.
Мы хотим описать временную эволюцию электронной структуры, индуцированной внешним потенциалом. Для этого рассмотрим сначала конструкцию молекулярной динамики Эрнфеста.
Так как мы можем рассчитать электронные орбитали (не только основное состояние, но и возбужденные состояния), то мы можем построить редуцированную модель, основанную на молекулярной
динамике Эрнфеста и временной теории функционала плотности [9]. Цель — получить возмущение электронной плотности из-за внешнего электрического потенциала и соответствующую силу, действующую на ядра водорода в центральной полости молекулы.
Разложим электронную волновую функцию по полному набору функций:
Ф(г,11,г) = ^2ст(1)фт(г,П),
т
где г и К — электронные и ядерные положения. Пусть фт{г, К) удовлетворяют стационарному уравнению
Нфт(г, К) = Етфт(г, К).
Молекулярная динамика Эрнфеста трактует временную зависимость коэффициентов ст(1) как определяющую эволюцию системы во времени. Запишем гамильтониан в следующей форме:
Н = Но + %н + Уех\
где Уен — потенциал электростатического взаимодействия между электронами и двумя ядрами водорода в центральной области молекулы, На — полный молекулярный гамильтониан с исключенным V, а. К"*' — внешний потенциал.
Мы вводим аппроксимацию невозмущенной потенциальной поверхности молекулы следующим образом:
Up = -bpeq; UF(s) = cos ( —)
4 \ iimoi J
где peq — "равновесная" электронная плотность полной молекулы. Ее профиль вдоль координационной переменной вычислен нами с помощью кода CPMD.
Рассмотрим движение только двух ядер водорода (описываемое одной координационной переменной), возбужденное электронной плотностью из-за действия приложенного внешнего потенциала. Мы вводим фиктивную частицу, которая показывает положение ядер водорода как s(t), и получим приближенные уравнения движения в следующей форме:
= ~ Iе™!2 "ЕЕ cm(t)ck(t)(Em - Ек) dmk(s), (2)
т>гп0 т>гп0 кфт,
= Cm{t)Em{s) ■ dmk(s) + (з)
где
dm^k(s(t)) = (ф,
ds
Фк
Вероятность конечного состояния системы дает |Ф|2. Рассмотрим орбитали HOMO, LUMO и LUMO + 1, включенные в разложение волновой функции. Структура орбиталей LUMO и LUMO + 1 показана на рис. 1, в, г.
Мы выбирали внешний потенциал в следующей форме:
Fext = VQ ■ (z - zo)e{r-ro)2/^ cos(wi).
Решение системы уравнений (2) и (3) показано на рис. 3, а ,б. Здесь представлена временная зависимость коэффициентов cm(t) для двух значений внешнего приложенного потенциала: Vq = 0.2 eV (рис. 3, а) и Vq = 0.5 eV (рис. 3,6); для НОМО-орбитали — сплошные линии, LUMO-орбитали — точечная линия, (LUMO + 1)-орбитали — штриховая линия. Видно, что величина Vq = 0.2 eV недостаточна для переключения, и переключение существует, когда Vq = 0.5 eV. Мы можем наблюдать частоту переключения, используя различную амплитуду и частоту внешнего потенциала, потому что процесс может рассматриваться как переключаемый благодаря своей пространственной структуре.
Важно заметить, что переключение может наблюдаться, только если имеется резонанс между потенциалом и частотой перехода между орбиталями и разделение проводится амплитудой Уо сигнала. Критическое значение приложенного потенциала, необходимого для переключения, соответствует найденному методом метадинамики.
Рис. 3. Временная зависимость коэффициентов для двух значений приложенного потенциала: Vo = 0.2 eV (а) и Vo = 0.5 eV (б); для НОМО-орбитали — сплошные линии, LUMO-орбитали — точечная линия, (LUMO + 1)-орбитали — штриховая линия
5. Заключение. Представлено моделирование "из первых принципов" молекулярного переключения на основе индуцированной током реакции изомеризации молекулы нафталоцианина. С помощью численного кода CPMD рассчитан профиль поверхности свободной энергии и путь реакции. Вычис-
ления проведены на супер-ЭВМ IBM Blue Gene/P. Найдена высота энергетического барьера реакции вдоль координационного направления. Показано, что динамика всех ядер в молекуле (колебательные и вращательные моды) важна для нахождения реакционного пути.
На основе вычислений "из первых принципов" предложена редуцированная модель Эрнфеста, позволяющая рассматривать самосогласованную эволюцию возбужденных электронных состояний и динамики ядер водорода. Найдено, что величина потенциала, необходимая для переключения, находится в хорошем соответствии с высотой барьера поверхности свободной энергии, полученной методом ме-тадинамики.
СПИСОК ЛИТЕРАТУРЫ
1. Rieth M.,Schommers W. Handbook of theoretical and computational nanotechnology. Germany: Forschungszentrum Karlsruhe, 2006.
2. Попов A.M. Вычислительные нанотехнологии. M.: МАКС Пресс, 2009.
3. Liljeroth P., Repp J., Meyer G. Current-induced hydrogen tautomerization and conductance switching of naphthalocyanine molecules // Science. 2007. 317. P. 1203-1206.
4. Car R., Parrinello M. Unified approach for molecular dynamics and density-functional theory // Phys. Rev. Lett. 1985. 55. P. 2471-2474.
5. Andreoni W., Curioni A. New advances in chemistry and materials science with CPMD and parallel computing // Parallel Computing. 2000. 26. P. 819-842.
6. Ensing В., Laio A., Parrinello M. et al. A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions //J. Phys. Chem. B. 2005. 109. N 14. P. 6676-6687.
7. Kohn W. Density functional and density matrix method scaling linearly with the number of atoms // Phys. Rev. Lett. 1996. 76. N 17. P. 3168-3171.
8. Fukui K. The path of chemical reactions. The IRC approach // Acc. Chem. Res. 1981. 14. P. 363-368.
9. Andrade X., Castro A., Zueco D. et al. A modified Ehrenfest formalism for efficient large-scale ab initio molecular dynamics //J. Chem. Theory Comput. 2009. 5. N 4. P. 728-742.
Поступила в редакцию 14.01.10
MULTISCALE MODEL OF MOLECULAR SWITCH
Popov A. M.. Shumkin G. N.
Ab-initio simulations of molecular switch on the base of STM induced isomerization of naphthalocyanine molecule is presented. Free energy profile and reaction path were analyzed using metadynamics method in the frame of Car-Parrinello Molecular Dynamics (CPMD) code. Simulations were done on IBM Blue Gene/P supercomputer at Moscow State University. Barrier height for reaction along coordination variable is found for isomerization reaction. Multiscale model of molecular switch is proposed on the base of Ehrenfest Molecular Dynamics and the excited electronic states. Model uses coordination direction of reaction and free energy profile, obtained in the detailed CPMD calculations of whole molecule structure.
Keywords: molecular switches, quantum molecular dynamics, computational multiscale model, chemical reaction path.
