Файл: Исследование свойств углеродных нанотрубок методами молекулярной динамики.pdf
ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 29.10.2023
Просмотров: 18
Скачиваний: 2
ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.
Исследование свойств углеродных нанотрубок
методами молекулярной динамики
Д.В. Зленко a,1
, П.А. Мамонов a,1
, А.М. Нестеренко a,2
a
Группа теоретической биофизики: http://erg.biophys.msu.ru/
1
Кафедра биофизики Биологического факультета МГУ им. М.В. Ломоносова, Москва, Россия
2
НИИ Физико-химической биологии им. А.Н. Белозерского МГУ им. М.В. Ломоносова, Москва, Россия
Аннотация
В рамках данной задачи студенты получают практические навыки создания и настройки молекулярно-механических моделей, осваивают методы молекулярного моделирования. В качестве объекта моделирования используются углеродные нанотрубки различного строения. В ходе работы студенты исследуют их механические колебательные свойства, упругие свойства, а также проницаемость нанотрубок для молекул воды. При решении этих задач студенты осваивают методы гармонического анализа, управляемой и обычной молекулярной динамики, знакомятся с техникой интерпретации результатов молекулярного моделирования. Перечисленный комплекс работ позволяет дать студентам минимальный базовый уровень владения методами молекулярного моделирования, а также познакомить с интересными и необычными свойствами углеродных нанотрубок.
Введение.
На сегодняшний день методы молекулярного моделирования стали незаменимым инструментом исследования биологических молекул и наноструктур. Совершенствование теоретических подходов и развитие вычислительной техники позволило им приблизиться по информативности к экспериментальным методам исследования. Вместе с тем подходы
1
основанные на моделировании обладают рядом преимуществ перед экспериментальными методами. Исследование свойств молекулярных систем «in silico» не требует специализированного экспериментального оборудования и позволяет в короткие сроки решать задачи, требующие значительного времени и средств при исследовании в реальном эксперименте. Несмотря на то, что в общем случае точность результатов, получаемых при моделировании, уступает экспериментальным, подобный подход находит применение в ситуации, когда стоит задача массового тестирования свойств большого количества соединений. Например, таким образом тестируется сродство потенциальных лигандов к белку-мишени, при создании новых лекарств. Другим преимуществом подхода основанного на моделировании является возможность более глубокого анализа процессов происходящих в молекулярной системе. Более того, в таких сложных системах, как биополимеры и их комплексы, а также наноструктуры, интерпретация результатов реальных экспериментов часто оказывается затруднительной без привлечения результатов численного моделирования.
Одной из таких задач является, например, исследование структурно-функциональных взаимосвязей в белках, и в частности механизмов конформационной регуляции их активности. Наконец, подходы основанные на моделировании позволяют исследовать свойства еще не созданных, гипотетических био- и наноструктур, что может привести к созданию новых материалов и наноустройств.
Одним из широко применяемых подходов к моделированию молекулярных систем, знакомству с которым посвящен данный практикум, является молекулярно-механический подход. В основе этого подхода лежит представление об атомах, как о классических точечных взаимодействующих частицах. Взаимодействие между этими частицами описывается набором эмпирических потенциалов, именуемых «силовым полем». Наиболее часто применяемый для моделированию молекулярных систем метод состоит в расчете динамики молекулярно-механической системы в рамках классической механики, с последующим анализом полученных траекторий движения атомов. Такой подход обладает высокой информативностью, позволяя описывать большинство процессов, происходящих в молекулярной системе и не сопровождающихся разрывом/образованием ковалентных связей и изменением электронного состояния системы. Важным преимуществом молекулярно- механического моделирования является относительно низкая ресурсоемкость соответствующих вычислительных процедур, что позволяет рассчитывать динамику систем включающих сотни тысяч атомов на временах до десятков наносекунд.
В качестве объекта исследования в практикуме используются нанотрубки различного строения. Нанотрубки являются классическим нанообъектом и представляют интерес как с
2
теоретической, так и с практической точек зрения. Данный практикум посвящен исследованию механических свойств нанотрубок и их проницаемости для молекул воды.
Первой задачей, решаемой в рамках практикума, является исследование спектра нормальных колебаний изолированной нанотрубки. Данная задача призвана продемонстрировать то промежуточное положение между макрообъектами и молекулярными системами, которое занимают нанообъекты. Так, будучи одиночной молекулой, нанотрубка демонстрирует свойства присущие механическим макросистемам: в колебательном спектре отдельной нанотрубки присутствуют колебания, аналогичные колебаниям струны, деформации сжатия-растяжения и изгибной деформации упругой балки.
Следующая задача практикума состоит в определении модуля упругости сжатия отдельной нанотрубки. С этой целью выполняется расчет управляемой молекулярной динамики нанотрубки при действии на нее растягивающей силы до достижения системой состояния равновесия. По изменению длины нанотрубки и значению приложенной силы рассчитывается модуль Юнга. Также рассчитывается модуль упругости материала на основе нанотрубок с учетом плотности их упаковки. Результаты полученные на данном этапе демонстрируют более высокие прочностные характеристики нанотрубок по сравнению с известными конструкционными материалами.
Заключительная часть практикума посвящена моделированию латеральной проводимости нанотрубок для воды. С этой целью выполняется расчет молекулярной динамики нанотрубок в воде и вычисляются коэффициенты диффузии молекул воды в полости нанотрубки и в объеме воды. Результаты моделирования свидетельствуют о высокой проводимости нанотрубок для воды. Способность нанотрубок к латеральной проводимости веществ считается практически важным свойством, которое может быть использовано для организации трансмембранной проницаемости липидных мембран, при решении ряда медико-биологических проблем.
Цель и задачи
Основной целью данного практикума является обучение студентов принципам создания молекулярно-механических моделей, подбору параметров проведения молекулярных расчетов, а также основным приемам моделирования молекулярных систем и обработки полученных массивов данных.
В ходе выполнения практикума студенты решают следующие задачи:
1. Создание молекулярно-механических моделей углеродных нанотрубок различного
3
строения.
2. Оптимизация геометрии молекулярных моделей нанотрубок.
3. Расчет и анализ спектра нормальных колебаний нанотрубок.
4. Расчет модуля упругости нанотрубки при растяжении, с использванием метода управляемой молекулярной динамики.
5. Расчет коэффициента диффузии воды внутри нанотрубки и в растворе с использованием метода молекулярной динамики.
Объекты и методы
В качестве объекта исследований в задаче используются углеродные нанотрубки различного строения. Моделирование осуществляется в рамках молекулярно-механического подхода. Для моделирования используется свободно распространяемый пакет программ
GROMACS [http://www.gromacs.org/]. Для создания, редактирования и визуализации молекулярных систем используется программа PyMOL. Непосредственно для создания нанотрубок и визуализации нормальных колебаний используются программные модули
PyMOL ntgen.py и nmgen.py, доступные на сайте http://erg.biophys.msu.ru/. В качестве среды для технических расчетов используется интерактивная командная среда IPython
[http://ipython.scipy.org/], являющаяся интерфейсом интерпретатора языка Python
[http://www.python.org/]. В сочетании с библиотеками численных алгоритмов SciPy
[http://www.scipy.org/] и графики matplotlib [http://matplotlib.sourceforge.net/] данное ПО является мощным и удобным пакетом для численного анализа.
Ход задачи практикума.
Часть 1. Создание нанотрубок
В качестве объекта исследования в задаче предлагается использовать нанотрубки различного строения. Нанотрубку можно представить как прямоугольный фрагмент монослоя графита, свернутый в цилиндр (рис. 1). Строение нанотрубки определятся ориентацией сторон этого прямоугольника относительно векторов трансляционной симметрии монослоя графита (a
1
, a
2
). Практически, строение нанотрубки задается парой индексов (n,m), определяющих одну из сторон этого прямоугольника (C
h
) как линейную комбинацию векторов a
1
и a
2
: C
h
= n a
1
+ m a
2
, при этом вектор перпендикулярный C
h оказывается параллелен оси получающегося при сворачивании цилиндра. В данной задаче предлагается использовать две нанотрубки: строения (10,0), условно называемая «зигзаг», и
4
(6,6) - «кресло» (рис. 2). Данный выбор позволяет продемонстрировать существенное влияния строения нанотрубок на их механические и электрические свойства.
Практически, генерация координат нанотрубок осуществляется при помощи специального программного модуля PyMOL. Для генерации координат следует запустить программу PyMOL и выполнить из ее командной строки следующие команды:
PyMOL >
ntgen nt100, 10, 0, 24, save=nt100
PyMOL >
ntgen nt66, 6, 6, 20, save=nt66
В результате будут сгенерированы объекты nt-10-0 и nt-6-6 с нанотрубками типа
«зигзаг» и «кресло» соответственно (рис. 2). Также в текущей директории будут созданы соответствующие файлы со структурой нанотрубок (*.gro) и с их топологией (*.top) пригодные для дальнейшего использования с пакетом GROMACS.
Дальнейшие манипуляции выполняются для каждой из полученных структур, а при их описании будет использовано единое обозначение nt-N-M, подразумевающее обе структуры.
Часть 2. Оптимизация геометрии системы
5
Рис. 1: Развертка нанотрубки на монослое графита. a1, a2 — вектора
трансляционной симметрии монослоя графита. Ch — сторона фрагмента
монослоя графита, изгибающаяся при сворачивании цилиндра. Задана как
линейная комбинация векторов трансляционной симметрии a1, a2 с
коэффициентами n и m. T — вектор направленный вдоль оси симметрии
нанотрубки, и задающий ее длину.
Чтобы перейти к расчетам молекулярной динамики, нам необходимо создать основной файл топологии (*.top). Он должен ссылаться на файлы молекулярных топологий (*.itp) и подключать таблицы силового поля (ff*.itp). Также файл должен содержать информацию о том, в каком количестве представлены в структуре молекулы и в каком порядке они идут.
Файл составляется в текстовом редакторе и выглядит следующим образом:
#include "ffoplsaa.itp"
#include "nanotube.itp"
[ system ]
Your nanotube name
[ molecules ]
CNT 1
В последней строчке указывается, что в нашей системе имеется один остаток CNT.
Следующий этап - создание файлов молекулярно-динамических параметров (*.mdp), в которых содержаться все условия проведения МД-расчетов: один - для оптимизации геометрии (em.mdp), а второй - для расчета матрицы Гесса исследуемой нанотрубки (nm.mdp
). Эти файлы студенты создают вручную, под руководством сотрудника проводящего практикум. Ниже дан пример файла em.mdp с краткими комментариями, предназначенного для проведения оптимизации геометрии.
;
em.mdp
cpp
=
/usr/bin/cpp ;
Препроцессор.
integrator
=
steep
;
Метод оптимизации геометрии. nsteps
=
100000
;
Кол-во шагов интегрирования.
emtol
=
100
;
Порог энергии.
6
Рис. 2: Нанотрубки различного строения. Водороды на конца нанотрубок не
показаны. а) Нанотрубка строения (10,0), «зигзаг». б) Нанотрубка строения
(6,6), «кресло».
emstep
=
0.001
;
Величина шага по энергии.
nstcomm
=
1
;
Частота удаления движения ЦМ.
comm_mode =
linear
;
Способ удаления движений ЦМ.
ns_type
=
grid ;
Способ поиска соседей.
nstlist
=
1
;
Частота обновления списка соседей rlist
=
1.25
;
Радиус обрезания списка соседей.
coulombtype
=
pme ;
Способ расчета электростатики.
fourierspacing =
0.12 ;
Амплитуда функций в Фурье преобр.
pme_order =
4
;
Порядок интерполяции.
ewald_rtol
=
1e-5
;
Точность расчета электростатики.
ewald_geometry =
3d
;
Размерноть пространства.
optimize_fft
=
yes ;
Оптимизация fft.
rcoulomb =
1.25 ;
Радиус обрезания электростатики. vdwtype
=
switch
;
Способ описания VdW взаимодействий.
rvdw =
1.15 ;
Радиус обрезания VdW
rvdw_switch
=
1.0
;
Радиус переключения режимов VdW
7
Рис. 3: Нанотрубки с различными типами симметрии.
Вверху - "кресло", посередине - "зигзаг", внизу -
хиральная нанотрубка. Хорошо видны различия в
расположении колец и строении края нанотрубки.
pbc
=
xyz
;
Периодичесие граничные условия.
nstxout
=
10
;
Частота вывода координат.
nstlog
=
10
;
Частота вывода записей в лог-файл.
nstenergy =
10
;
Частота вывода энергий.
Затем проводится оптимизация геометрии построенных на первом этапе структур в вакууме. Все расчеты GROMACS запускаются в 2 этапа: сбор всей информации по расчету в один бинарный файл (GROMacsPreProcessor, grompp) и собственно численное интегрирование уравнений движения (MolecularDynamicsRun, mdrun). На первом этапе создается полностью самостоятельный бинарный файл, который необходим и достаточен для запуска расчета filename.tpr. Программа mdrun работает, получая на вход только этот файл.
Запуск команды производится следующим образом:
$ grompp f em.mdp c *.gro p *.top o em.tpr
Программа mdrun дает на выходе от 4 различных файлов: лог расчета (mdrun.log), бинарный файл макроскопических характеристик системы (ener.edr), бинарный файл координат, скоростей и сил всех атомов системы (traj.trr), файл с состоянием структуры на последний момент расчета (confout.gro). Для того, чтобы не задавать имена каждому файлу используется директива -deffnm.
$ mdrun s em.tpr deffnm common_filename
Следует обратить внимание на тот факт, что конечная структура confout.gro имеет исключительно декоративный смысл, потому как не может быть использована в дальнейших расчетах из-за низкой точности. Конечное состояние системы с наивысшей точностью записано только внутри файла traj.trr и именно он должен использоваться для дальнейших запусков. Это имеет большое значение для корректного расчета Гессиана системы, необходимого для анализа нормальных колебаний системы.
Часть 3. Расчет нормальных колебаний
Для расчета матрицы Гесса также используется утилита mdrun, на вход ей необходимо подать траекторию, полученную в ходе оптимизации:
$ grompp f nm.mdp c *.gro t enmin.trr p *.top o nm.tpr
Файл nm.mdp должен отличаться от предыдущего файла параметров только значением integrator = nm.
$ mdrun s nm.tpr deffnm hessename
Полученный на выходе Гессиан (файл hessename.mtx) представляет из себя матрицу гессе в бинарном виде, ее необходимо передать в программу для pyMol — nmgen
8
(
http://erg.biophys.msu.ru/piton/nmgen.py
). Ее таже необходимо подключить командой import, как это делалось в случае ntgen. Программа рассчитывает частоты и вектора нормальных колебаний, а также позволяет визуализировать интересующие нормальные моды и непосредственно наблюдать за смещениями вовлеченных в эти колебания атомов. Пример такой визуализации представлен на Рис.2.
> nmgen name, id, amp=n где name - имя нанотрубки, для которой рассчитывается визуализация колебания, id - номер этого колебания, а n - амплитуда смещения атомов. Необходимо отыскать среди низкочастотных нормальных колебаний, колебания соответствующие деформации нанотрубки как целого, выписать их частоты и описать эти колебания. Следует обратить внимание на то, как меняются частоты нормальных колебаний при изменении числа вовлеченных в них атомов, а также при изменении кратности (числа минимумов и максимумов) однотипных колебаний, например, деформации сжатия-растяжения нанотрубки.
Часть 4. Расчет модуля Юнга углеродной нанотрубки
Существут несколько способов вычисления модуля Юнга из МД-расчета. Первые два способа основаны на результатах управляемой молекулярной динамики, при котором к сумме всех сил, действующих на краевые атомы нанотрубки добавляется некоторая внешняя сила, при помощи которой моделируется ее растяжение (Рис.4). Для вычисления модуля Юнга используем уравнение (Формула 1).
Рассчитать модуль Юнга необходимо двумя способами, что оправданно, так как часть энергии, привнесенной внешней силой, расходуется на разогрев системы и теряется при термостатировании. Это можно учесть, если рассматривать непосредственно потенциальную энергию. Соответственно, первый способ вычисления дает завышенный результат, что наглядно демонстрирует особенности метода управляемой молекулярной динамики.
9
E = (F/S)/(
∆
l/l)
Формула 1: Здесь F — внешняя сила, S — площадь сечения нанотрубки,
∆
l —
удлинение, l — начальная длина. Можно задать силу F явно в mdp-файле и
подставить это значение в формулу или же подставить в качестве силы
F=
∆
E/
∆
l, где
∆
E — разница потенциальной энергии системы до и после
расяжения.
Рис. 4: Растяжение нанотрубки под действием
внешней силы. Расчеты проводились методом
управляемой молекулярной динамики.
Расчет модуля Юнга с помощью МД под действием внешней силы происходит по следующей схеме:
1.
Термостатирование системы в отсутствии внешней силы (200 пс).
2.
Термостатирование системы с приложенной внешней силой (200 пс).
3.
Визуальная проверка расчетов на корректность.
4.
Получение средних значений удлинения и энергии.
Необходимо создать два файла молекулярно-динамических параметров, один из которых предназначен для расчетов управляемой молекулярной динамики (force.mdp), а
10
Рис. 5: Примеры нормальных колебаний нанотрубки.
Для наглядности амплитуды сильно преувеличены.
второй - для проведения термостатирования (md.mdp). Файл параметров для расчетов динамики существенно отличается от предыдущих, поэтому рассмотрим его здесь подробнее.
cpp = /usr/bin/cpp integrator = sd — интегратор молекулярной динамики dt = 0.0007
— шаг интегрирования ур-й движения (пс)
nsteps = 100000000 — кол-во шагов nstcomm = 10
comm_mode = Linear ld_seed = 12173 — иницализатор случайных чисел nstxout = 2000 — частота записи координат nstlog = 500 — частота записи в log-файл nstenergy = 5000 — частота записи в бинарный файл энергий energygrps = system nstlist = 10
ns_type = grid
; это параметры вычисления энергии — скопируйте их из предыдущего mdp-файла pbc = *
rlist = *
coulombtype = *
fourierspacing = *
pme_order = *
ewald_rtol = *
ewald_geometry = * optimize_fft = *
rcoulomb = *
vdwtype = *
rvdw = *
rvdw_switch = *
; настройки термостата и баростата tcoupl = berendsen ; тип баростата tc_grps = system tau_t = 0.1 ; частота срабатывания баростата ref_t = 280 ; температура pcoupl = berendsen pcoupltype = no ; gen_vel = yes gen_temp = 280 ; начальная температура gen_seed = 15000
; ускорение (только для force.mdp)
acc_grps
= ltail htail accelerate
= 0 0 -5 0 0 +5
Запуск расчетов происходит по той же схеме:
$ grompp -f force.mdp -c *.gro -p *.top -o force.tpr -n *.ndx
11
$ mdrun -s force.tpr -deffnm force
$ grompp -f md.mdp -c *.gro -p *.top -o md.tpr
$ mdrun -s md.tpr -deffnm md
Обратите внимание на то, что для подготовки одного из tpr-файлов необходим файл
NDX. Как его составить, смотрите ниже. На выходе получаются траектория (*.trr) - файл, в котором содержаться координаты всех атомов системы во все моменты времени и энергетический файл (*.edr).
По окончании молекулярно-динамических расчетов можно приступать к обработке полученных данных и расчету модуля Юнга. Для этого на первом этапе необходимо рассчитать изменение длины и потенциальной энергии нанотрубки в релаксированном и растянутом состоянии.
Рис. 6: Динамика изменения потенциальной энергии
нанотрубки при приложении к ней растягивающей силы.
Для этого используются соответствующие утилиты пакета GROMACS: g_traj и g_energy. Синтаксис g_traj следующий:
$ g_traj -f *.trr -s *.tpr -n *.ndx -ox coord.xvg -com -ng 2 -nox
-noy
[-com] означает, что будут выведениы координаты центра масс группы, [-ng 2] означает,
12
что будут обрабатываться 2 группы атомов, а ключи [-nox/-noy] означают, что координаты X и
Y выводится не будут. Для того, чтобы указать утилите g_traj координаты каких именно атомов необходимо извлечь из траектории используется индекс-файл (*.ndx), который создается программой make_ndx и далее редактируется вручную с использованием пакета
PyMOL. Для генерации начального индекса, запустите следующую команду и нажмите q-
Enter:
make_ndx -f file.gro -o indexname.ndx
В полученный NDX-файл необходимо добавить две группы атомов: верхнее углеродное кольцо и нижнее. Синтаксис следующий:
[htail]
1 2 3 4 5 6 7 8
[ltail]
467 456 778 346 758
Здесь htail и ltail — название групп, а номера — это абсолютные номера атомов, входящих в группу. htail — группа концевых углеродов лежащая «вверху» по оси Z, ltail — то же, но с другой стороны нанотрубки.
Синтаксис команды g_energy следующий:
$ g_energy -f *.edr -s *.edr -o ener.edr
При просьбе выбрать поля из таблице, необходимо выбрать полную потенциальную и кинтечискую энергии.
На выходе обоих программ получаются файлы XVG - текстовые файлы с таблицами интересующих значений, которые затем обрабатываются в среде IPython, позволяющей, в частности визуализировать изменения наблюдаемых параметров во времени (Рис.5).
Следующий этап работы: вычисление модуля Юнга на основе имеющихся таблиц энергий и координат. Вычисления производятся в MATLAB-подобной среде IPython. Для загрузки XVG-файла в массив используется функция load:
> a = load('file.xvg')
Для вычисления среднего значения массива — функция mean:
> eav = mean(a[:,3])
Для визуализации массива — функция plot:
13
> plot(a[:,0], a[:,1])
Вычисления модуля Юнга производятся согласно формуле 1. В одном случае в качестве силы F подставляется масса одного углеродного кольца, умноженного на ускорение, заданное в mdp-файле. Во втором случае в качестве F подставляется разница средних значений потенциальной энергии системы, деленная на среднее удлинение по оси Z.
Обратите внимание, что в обоих случаях необходимо брать среднее значение энергии и удлинения для равновесной области. Итерация, начиная с которой систему можно считать равновесной определяется качественно по графику энергии (см. Рис. 6).
Для вычислений необходимо учитывать размерности физических единиц, принятые в
GROMACS:
1. Энергия — кДж/моль
2. Ускорение — нм/пс
2 3. Расстояние — нм
4. Время — пс
Часть 3. Скорость диффузии воды в углеродной нанотрубке
Рассчитать коэффициент диффузии, исходя из достаточно коротких (200-500 пс) траекторий корректно невозможно из-за существенного вклада в суммарное смещение центров масс каждой конкретной молекулы случайных тепловых движений молекул.
Поэтому необходимо нормировать полученное для коэффициента диффузии воды внутри нанотрубки (Рис.7,8) значение, используя коэффициент диффузии воды в объемной фазе, так как этот параметр, при использовании более длинных (10-20 нс) траекторий хорошо согласуется с известным из эксперимента значением.
Для выполнения этой части работы необходимо сначала добавить в систему воду. Эту операция выполняется при помощи утилиты gen_box:
$ gen_box -sp *.gro -cs tip4p.gro -o *.gro
Как нетрудно видеть, для модели воды используется модель tip4p, в которую добавлена четвертая виртуальная частица, описывающая неподеленную электронную пару атома кислорода, что необходимо для корректного образования водородных связей в системе.
Файл молекулярно-динамических параметров подготавливается по уже известной схеме. Файлы TOP необходимо модифицировать в связи с добавлением в систему воды, а
14
именно добавить строчку ..
SOL 256
.. в TOP-файл. После окончания молекулярно динамических расчетов необходимо извлечь из полученной траектории координаты молекул воды, которые необходимы для расчета коэффициента диффузии. Для этого используется утилита g_traj. Соответственно, ей необходимо передать номера атомов интересующих молекул воды, что делается при помощи соответствующего индекс-файла (*.ndx). Создание индекс-файла уже обсуждалось ранее. В данном случае необходимо внести группу молекул воды, находящихся внутри нанотрубки.
Затем из траектории извлекаются координаты атомов интересующих молекул воды:
$ g_traj -f *.trr -s *.tpr -n *.ndx -ox coord.xvg -mol
, где ключ [-mol] указывает на то, что будут выписаны не координаты атомов, а координаты центров масс молекул.
Обработка полученных данных производится также в среде IPython. Написать программу необходимо самостоятельно, руководствуясь формулой:
D = (X(t)-X(0))^2/2nt
.., где D – коэффициент диффузии, X(t) — радиус-вектор положения молекулы в момент t, n – размерность пространства. Параметр n будет отличаться для молекул воды в объемной фазе и внутри нанотрубки.
Требования к отчету
В конце практикума студенты готовят письменный отчет по задаче. Отчет должен включать следующие разделы:
−
Введение, включающее краткую характеристику метода молекулярной динамики, а также характеристика объекта исследования - нанотрубок, их структуры, свойств и возможных путей их практического использования.
−
Цель и основные задачи данного практикума.
−
Описание хода выполнения задачи, состоящее их трех частей, посвященных трем разделам задачи. Каждая часть должна начинаться с постановки задачи и краткого описания используемых для ее решения подходов. Затем должно следовать собственно описание хода работы, снабженное промежуточными числовыми
15
результатами и иллюстрациями. В заключении приводятся финальные результаты моделирования, со всеми необходимыми выкладками, сравниваются и интепретируются результаты моделирования для нанотрубок различного строения.
−
Заключение, в котором кратко формулируются полученные результаты и дается сравнительная характеристика нанотрубок различного строения.
Контрольные вопросы
1. Сколько атомов содержит минимальная кристаллографическая ячейка монослоя графита? Пояснить на схеме.
2. Каковы типичные размеры (радиус, длина) нанотрубки?
3. Рассчитать погрешности в определении равновесной длины связи в двухатомной молекуле и ее потенциальной энергии, при оптимизации с пороговым значением градиента энергии 10 кДж моль
-1
нм
-1
. Потенциальная энергия системы задана, как
U = k (R
AB
-R
0
)
2
/ 2.
4. Справедливы ли результаты гармонического анализа для системы не находящейся в потенциальном минимуме поверхности потенциальной энергии? Пояснить на примере одномерного потенциала.
5. Найти частоты и вектора нормальных колебаний одномерной системы состоящей из двух частиц массой m
1
и m
2
, если потенциальная энергия данной системы задана как
U = k
1
(x
1
-x
2
)
2
/ 2.
6. Что такое шаг интегрирования и радиус обрезания в методе молекулярной динамики?
16
Рис. 7: Цепочки связанных друг с другом водородными связями
молекул воды внутри нанотрубки.
Чем важны эти параметры?
7. Что такое периодические граничные условия в методе молекулярной динамики, для чего они используются?
8. В чем состоит метод управляемой молекулярной динамики?
9. Что такое модуль упругости? Какие существуют методы расчета коэффициента упругости на основании данных молекулярной динамики?
10. Что такое коэффициент диффузии? Какие микроскопические характеристики молекулярной системы его определяют?
11. Почему коэффициент диффузии воды внутри нанотрубки выше, чем коэффициент диффузии воды в объемной фазе?
Рекомендуемая литература
1. Сивухин Д.В. Общая физика. Москва. Наука. 1974.
2. Тихонов А.Н., Самарский А.А. Уравнения математической физики. Москва. Наука.
1972.
3. David van der Spoel, Lindahl E., Hess B., el. all. GROMACS User Manual. Department of
Biophysical Chemistry, University of Groningen. Nijenborgh 4, 9747 AG Groningen, The
17
Рис. 8: Молекулы воды внутри нанотрубки. Хорошо видно, что
молекулы воды располагаются внутри нанотрубки в весьма
ограниченном объеме.
Netherlands. 1991-2000.
4. Suzuki K. On Elastic Properties of Single-walled Carbon Nanotubes as Composite
Reinforcing Fillers. J. Comp. Mat. 41(9):1123-1135. 2007.
5. ChenXin F.U., YunFei C., JiWei J. Molecular dynamics simulation of the test of single-
walled carbon nanotubes under tensile loading. Sci China Ser E-Tech Sci. 50(1):1-17. 2007.
6. Garde A.S., Hummer G. Osmotic water transport through carbon nanotube membranes.
PNAS, 100(18):10175-10180. 2003.
7. Walthery J.H., Jaez R., Haliciogluz T., Koumoutsakosy P. Molecular dynamics simulations
of carbon nanotubes in water. Center for Turbulence Research Proceedings of the Summer
Program 2000 18
Исследование свойств углеродных нанотрубок
методами молекулярной динамики
Д.В. Зленко a,1
, П.А. Мамонов a,1
, А.М. Нестеренко a,2
a
Группа теоретической биофизики: http://erg.biophys.msu.ru/
1
Кафедра биофизики Биологического факультета МГУ им. М.В. Ломоносова, Москва, Россия
2
НИИ Физико-химической биологии им. А.Н. Белозерского МГУ им. М.В. Ломоносова, Москва, Россия
Аннотация
В рамках данной задачи студенты получают практические навыки создания и настройки молекулярно-механических моделей, осваивают методы молекулярного моделирования. В качестве объекта моделирования используются углеродные нанотрубки различного строения. В ходе работы студенты исследуют их механические колебательные свойства, упругие свойства, а также проницаемость нанотрубок для молекул воды. При решении этих задач студенты осваивают методы гармонического анализа, управляемой и обычной молекулярной динамики, знакомятся с техникой интерпретации результатов молекулярного моделирования. Перечисленный комплекс работ позволяет дать студентам минимальный базовый уровень владения методами молекулярного моделирования, а также познакомить с интересными и необычными свойствами углеродных нанотрубок.
Введение.
На сегодняшний день методы молекулярного моделирования стали незаменимым инструментом исследования биологических молекул и наноструктур. Совершенствование теоретических подходов и развитие вычислительной техники позволило им приблизиться по информативности к экспериментальным методам исследования. Вместе с тем подходы
1
Одной из таких задач является, например, исследование структурно-функциональных взаимосвязей в белках, и в частности механизмов конформационной регуляции их активности. Наконец, подходы основанные на моделировании позволяют исследовать свойства еще не созданных, гипотетических био- и наноструктур, что может привести к созданию новых материалов и наноустройств.
Одним из широко применяемых подходов к моделированию молекулярных систем, знакомству с которым посвящен данный практикум, является молекулярно-механический подход. В основе этого подхода лежит представление об атомах, как о классических точечных взаимодействующих частицах. Взаимодействие между этими частицами описывается набором эмпирических потенциалов, именуемых «силовым полем». Наиболее часто применяемый для моделированию молекулярных систем метод состоит в расчете динамики молекулярно-механической системы в рамках классической механики, с последующим анализом полученных траекторий движения атомов. Такой подход обладает высокой информативностью, позволяя описывать большинство процессов, происходящих в молекулярной системе и не сопровождающихся разрывом/образованием ковалентных связей и изменением электронного состояния системы. Важным преимуществом молекулярно- механического моделирования является относительно низкая ресурсоемкость соответствующих вычислительных процедур, что позволяет рассчитывать динамику систем включающих сотни тысяч атомов на временах до десятков наносекунд.
В качестве объекта исследования в практикуме используются нанотрубки различного строения. Нанотрубки являются классическим нанообъектом и представляют интерес как с
2
Первой задачей, решаемой в рамках практикума, является исследование спектра нормальных колебаний изолированной нанотрубки. Данная задача призвана продемонстрировать то промежуточное положение между макрообъектами и молекулярными системами, которое занимают нанообъекты. Так, будучи одиночной молекулой, нанотрубка демонстрирует свойства присущие механическим макросистемам: в колебательном спектре отдельной нанотрубки присутствуют колебания, аналогичные колебаниям струны, деформации сжатия-растяжения и изгибной деформации упругой балки.
Следующая задача практикума состоит в определении модуля упругости сжатия отдельной нанотрубки. С этой целью выполняется расчет управляемой молекулярной динамики нанотрубки при действии на нее растягивающей силы до достижения системой состояния равновесия. По изменению длины нанотрубки и значению приложенной силы рассчитывается модуль Юнга. Также рассчитывается модуль упругости материала на основе нанотрубок с учетом плотности их упаковки. Результаты полученные на данном этапе демонстрируют более высокие прочностные характеристики нанотрубок по сравнению с известными конструкционными материалами.
Заключительная часть практикума посвящена моделированию латеральной проводимости нанотрубок для воды. С этой целью выполняется расчет молекулярной динамики нанотрубок в воде и вычисляются коэффициенты диффузии молекул воды в полости нанотрубки и в объеме воды. Результаты моделирования свидетельствуют о высокой проводимости нанотрубок для воды. Способность нанотрубок к латеральной проводимости веществ считается практически важным свойством, которое может быть использовано для организации трансмембранной проницаемости липидных мембран, при решении ряда медико-биологических проблем.
Цель и задачи
Основной целью данного практикума является обучение студентов принципам создания молекулярно-механических моделей, подбору параметров проведения молекулярных расчетов, а также основным приемам моделирования молекулярных систем и обработки полученных массивов данных.
В ходе выполнения практикума студенты решают следующие задачи:
1. Создание молекулярно-механических моделей углеродных нанотрубок различного
3
2. Оптимизация геометрии молекулярных моделей нанотрубок.
3. Расчет и анализ спектра нормальных колебаний нанотрубок.
4. Расчет модуля упругости нанотрубки при растяжении, с использванием метода управляемой молекулярной динамики.
5. Расчет коэффициента диффузии воды внутри нанотрубки и в растворе с использованием метода молекулярной динамики.
Объекты и методы
В качестве объекта исследований в задаче используются углеродные нанотрубки различного строения. Моделирование осуществляется в рамках молекулярно-механического подхода. Для моделирования используется свободно распространяемый пакет программ
GROMACS [http://www.gromacs.org/]. Для создания, редактирования и визуализации молекулярных систем используется программа PyMOL. Непосредственно для создания нанотрубок и визуализации нормальных колебаний используются программные модули
PyMOL ntgen.py и nmgen.py, доступные на сайте http://erg.biophys.msu.ru/. В качестве среды для технических расчетов используется интерактивная командная среда IPython
[http://ipython.scipy.org/], являющаяся интерфейсом интерпретатора языка Python
[http://www.python.org/]. В сочетании с библиотеками численных алгоритмов SciPy
[http://www.scipy.org/] и графики matplotlib [http://matplotlib.sourceforge.net/] данное ПО является мощным и удобным пакетом для численного анализа.
Ход задачи практикума.
Часть 1. Создание нанотрубок
В качестве объекта исследования в задаче предлагается использовать нанотрубки различного строения. Нанотрубку можно представить как прямоугольный фрагмент монослоя графита, свернутый в цилиндр (рис. 1). Строение нанотрубки определятся ориентацией сторон этого прямоугольника относительно векторов трансляционной симметрии монослоя графита (a
1
, a
2
). Практически, строение нанотрубки задается парой индексов (n,m), определяющих одну из сторон этого прямоугольника (C
h
) как линейную комбинацию векторов a
1
и a
2
: C
h
= n a
1
+ m a
2
, при этом вектор перпендикулярный C
h оказывается параллелен оси получающегося при сворачивании цилиндра. В данной задаче предлагается использовать две нанотрубки: строения (10,0), условно называемая «зигзаг», и
4
(6,6) - «кресло» (рис. 2). Данный выбор позволяет продемонстрировать существенное влияния строения нанотрубок на их механические и электрические свойства.
Практически, генерация координат нанотрубок осуществляется при помощи специального программного модуля PyMOL. Для генерации координат следует запустить программу PyMOL и выполнить из ее командной строки следующие команды:
PyMOL >
ntgen nt100, 10, 0, 24, save=nt100
PyMOL >
ntgen nt66, 6, 6, 20, save=nt66
В результате будут сгенерированы объекты nt-10-0 и nt-6-6 с нанотрубками типа
«зигзаг» и «кресло» соответственно (рис. 2). Также в текущей директории будут созданы соответствующие файлы со структурой нанотрубок (*.gro) и с их топологией (*.top) пригодные для дальнейшего использования с пакетом GROMACS.
Дальнейшие манипуляции выполняются для каждой из полученных структур, а при их описании будет использовано единое обозначение nt-N-M, подразумевающее обе структуры.
Часть 2. Оптимизация геометрии системы
5
Рис. 1: Развертка нанотрубки на монослое графита. a1, a2 — вектора
трансляционной симметрии монослоя графита. Ch — сторона фрагмента
монослоя графита, изгибающаяся при сворачивании цилиндра. Задана как
линейная комбинация векторов трансляционной симметрии a1, a2 с
коэффициентами n и m. T — вектор направленный вдоль оси симметрии
нанотрубки, и задающий ее длину.
Чтобы перейти к расчетам молекулярной динамики, нам необходимо создать основной файл топологии (*.top). Он должен ссылаться на файлы молекулярных топологий (*.itp) и подключать таблицы силового поля (ff*.itp). Также файл должен содержать информацию о том, в каком количестве представлены в структуре молекулы и в каком порядке они идут.
Файл составляется в текстовом редакторе и выглядит следующим образом:
#include "ffoplsaa.itp"
#include "nanotube.itp"
[ system ]
Your nanotube name
[ molecules ]
CNT 1
В последней строчке указывается, что в нашей системе имеется один остаток CNT.
Следующий этап - создание файлов молекулярно-динамических параметров (*.mdp), в которых содержаться все условия проведения МД-расчетов: один - для оптимизации геометрии (em.mdp), а второй - для расчета матрицы Гесса исследуемой нанотрубки (nm.mdp
). Эти файлы студенты создают вручную, под руководством сотрудника проводящего практикум. Ниже дан пример файла em.mdp с краткими комментариями, предназначенного для проведения оптимизации геометрии.
;
em.mdp
cpp
=
/usr/bin/cpp ;
Препроцессор.
integrator
=
steep
;
Метод оптимизации геометрии. nsteps
=
100000
;
Кол-во шагов интегрирования.
emtol
=
100
;
Порог энергии.
6
Рис. 2: Нанотрубки различного строения. Водороды на конца нанотрубок не
показаны. а) Нанотрубка строения (10,0), «зигзаг». б) Нанотрубка строения
(6,6), «кресло».
=
0.001
;
Величина шага по энергии.
nstcomm
=
1
;
Частота удаления движения ЦМ.
comm_mode =
linear
;
Способ удаления движений ЦМ.
ns_type
=
grid ;
Способ поиска соседей.
nstlist
=
1
;
Частота обновления списка соседей rlist
=
1.25
;
Радиус обрезания списка соседей.
coulombtype
=
pme ;
Способ расчета электростатики.
fourierspacing =
0.12 ;
Амплитуда функций в Фурье преобр.
pme_order =
4
;
Порядок интерполяции.
ewald_rtol
=
1e-5
;
Точность расчета электростатики.
ewald_geometry =
3d
;
Размерноть пространства.
optimize_fft
=
yes ;
Оптимизация fft.
rcoulomb =
1.25 ;
Радиус обрезания электростатики. vdwtype
=
switch
;
Способ описания VdW взаимодействий.
rvdw =
1.15 ;
Радиус обрезания VdW
rvdw_switch
=
1.0
;
Радиус переключения режимов VdW
7
Рис. 3: Нанотрубки с различными типами симметрии.
Вверху - "кресло", посередине - "зигзаг", внизу -
хиральная нанотрубка. Хорошо видны различия в
расположении колец и строении края нанотрубки.
=
xyz
;
Периодичесие граничные условия.
nstxout
=
10
;
Частота вывода координат.
nstlog
=
10
;
Частота вывода записей в лог-файл.
nstenergy =
10
;
Частота вывода энергий.
Затем проводится оптимизация геометрии построенных на первом этапе структур в вакууме. Все расчеты GROMACS запускаются в 2 этапа: сбор всей информации по расчету в один бинарный файл (GROMacsPreProcessor, grompp) и собственно численное интегрирование уравнений движения (MolecularDynamicsRun, mdrun). На первом этапе создается полностью самостоятельный бинарный файл, который необходим и достаточен для запуска расчета filename.tpr. Программа mdrun работает, получая на вход только этот файл.
Запуск команды производится следующим образом:
$ grompp f em.mdp c *.gro p *.top o em.tpr
Программа mdrun дает на выходе от 4 различных файлов: лог расчета (mdrun.log), бинарный файл макроскопических характеристик системы (ener.edr), бинарный файл координат, скоростей и сил всех атомов системы (traj.trr), файл с состоянием структуры на последний момент расчета (confout.gro). Для того, чтобы не задавать имена каждому файлу используется директива -deffnm.
$ mdrun s em.tpr deffnm common_filename
Следует обратить внимание на тот факт, что конечная структура confout.gro имеет исключительно декоративный смысл, потому как не может быть использована в дальнейших расчетах из-за низкой точности. Конечное состояние системы с наивысшей точностью записано только внутри файла traj.trr и именно он должен использоваться для дальнейших запусков. Это имеет большое значение для корректного расчета Гессиана системы, необходимого для анализа нормальных колебаний системы.
Часть 3. Расчет нормальных колебаний
Для расчета матрицы Гесса также используется утилита mdrun, на вход ей необходимо подать траекторию, полученную в ходе оптимизации:
$ grompp f nm.mdp c *.gro t enmin.trr p *.top o nm.tpr
Файл nm.mdp должен отличаться от предыдущего файла параметров только значением integrator = nm.
$ mdrun s nm.tpr deffnm hessename
Полученный на выходе Гессиан (файл hessename.mtx) представляет из себя матрицу гессе в бинарном виде, ее необходимо передать в программу для pyMol — nmgen
8
(
http://erg.biophys.msu.ru/piton/nmgen.py
). Ее таже необходимо подключить командой import, как это делалось в случае ntgen. Программа рассчитывает частоты и вектора нормальных колебаний, а также позволяет визуализировать интересующие нормальные моды и непосредственно наблюдать за смещениями вовлеченных в эти колебания атомов. Пример такой визуализации представлен на Рис.2.
> nmgen name, id, amp=n где name - имя нанотрубки, для которой рассчитывается визуализация колебания, id - номер этого колебания, а n - амплитуда смещения атомов. Необходимо отыскать среди низкочастотных нормальных колебаний, колебания соответствующие деформации нанотрубки как целого, выписать их частоты и описать эти колебания. Следует обратить внимание на то, как меняются частоты нормальных колебаний при изменении числа вовлеченных в них атомов, а также при изменении кратности (числа минимумов и максимумов) однотипных колебаний, например, деформации сжатия-растяжения нанотрубки.
Часть 4. Расчет модуля Юнга углеродной нанотрубки
Существут несколько способов вычисления модуля Юнга из МД-расчета. Первые два способа основаны на результатах управляемой молекулярной динамики, при котором к сумме всех сил, действующих на краевые атомы нанотрубки добавляется некоторая внешняя сила, при помощи которой моделируется ее растяжение (Рис.4). Для вычисления модуля Юнга используем уравнение (Формула 1).
Рассчитать модуль Юнга необходимо двумя способами, что оправданно, так как часть энергии, привнесенной внешней силой, расходуется на разогрев системы и теряется при термостатировании. Это можно учесть, если рассматривать непосредственно потенциальную энергию. Соответственно, первый способ вычисления дает завышенный результат, что наглядно демонстрирует особенности метода управляемой молекулярной динамики.
9
E = (F/S)/(
∆
l/l)
Формула 1: Здесь F — внешняя сила, S — площадь сечения нанотрубки,
∆
l —
удлинение, l — начальная длина. Можно задать силу F явно в mdp-файле и
подставить это значение в формулу или же подставить в качестве силы
F=
∆
E/
∆
l, где
∆
E — разница потенциальной энергии системы до и после
расяжения.
Рис. 4: Растяжение нанотрубки под действием
внешней силы. Расчеты проводились методом
управляемой молекулярной динамики.
Расчет модуля Юнга с помощью МД под действием внешней силы происходит по следующей схеме:
1.
Термостатирование системы в отсутствии внешней силы (200 пс).
2.
Термостатирование системы с приложенной внешней силой (200 пс).
3.
Визуальная проверка расчетов на корректность.
4.
Получение средних значений удлинения и энергии.
Необходимо создать два файла молекулярно-динамических параметров, один из которых предназначен для расчетов управляемой молекулярной динамики (force.mdp), а
10
Рис. 5: Примеры нормальных колебаний нанотрубки.
Для наглядности амплитуды сильно преувеличены.
cpp = /usr/bin/cpp integrator = sd — интегратор молекулярной динамики dt = 0.0007
— шаг интегрирования ур-й движения (пс)
nsteps = 100000000 — кол-во шагов nstcomm = 10
comm_mode = Linear ld_seed = 12173 — иницализатор случайных чисел nstxout = 2000 — частота записи координат nstlog = 500 — частота записи в log-файл nstenergy = 5000 — частота записи в бинарный файл энергий energygrps = system nstlist = 10
ns_type = grid
; это параметры вычисления энергии — скопируйте их из предыдущего mdp-файла pbc = *
rlist = *
coulombtype = *
fourierspacing = *
pme_order = *
ewald_rtol = *
ewald_geometry = * optimize_fft = *
rcoulomb = *
vdwtype = *
rvdw = *
rvdw_switch = *
; настройки термостата и баростата tcoupl = berendsen ; тип баростата tc_grps = system tau_t = 0.1 ; частота срабатывания баростата ref_t = 280 ; температура pcoupl = berendsen pcoupltype = no ; gen_vel = yes gen_temp = 280 ; начальная температура gen_seed = 15000
; ускорение (только для force.mdp)
acc_grps
= ltail htail accelerate
= 0 0 -5 0 0 +5
Запуск расчетов происходит по той же схеме:
$ grompp -f force.mdp -c *.gro -p *.top -o force.tpr -n *.ndx
11
$ mdrun -s force.tpr -deffnm force
$ grompp -f md.mdp -c *.gro -p *.top -o md.tpr
$ mdrun -s md.tpr -deffnm md
Обратите внимание на то, что для подготовки одного из tpr-файлов необходим файл
NDX. Как его составить, смотрите ниже. На выходе получаются траектория (*.trr) - файл, в котором содержаться координаты всех атомов системы во все моменты времени и энергетический файл (*.edr).
По окончании молекулярно-динамических расчетов можно приступать к обработке полученных данных и расчету модуля Юнга. Для этого на первом этапе необходимо рассчитать изменение длины и потенциальной энергии нанотрубки в релаксированном и растянутом состоянии.
Рис. 6: Динамика изменения потенциальной энергии
нанотрубки при приложении к ней растягивающей силы.
Для этого используются соответствующие утилиты пакета GROMACS: g_traj и g_energy. Синтаксис g_traj следующий:
$ g_traj -f *.trr -s *.tpr -n *.ndx -ox coord.xvg -com -ng 2 -nox
-noy
[-com] означает, что будут выведениы координаты центра масс группы, [-ng 2] означает,
12
Y выводится не будут. Для того, чтобы указать утилите g_traj координаты каких именно атомов необходимо извлечь из траектории используется индекс-файл (*.ndx), который создается программой make_ndx и далее редактируется вручную с использованием пакета
PyMOL. Для генерации начального индекса, запустите следующую команду и нажмите q-
Enter:
make_ndx -f file.gro -o indexname.ndx
В полученный NDX-файл необходимо добавить две группы атомов: верхнее углеродное кольцо и нижнее. Синтаксис следующий:
[htail]
1 2 3 4 5 6 7 8
[ltail]
467 456 778 346 758
Здесь htail и ltail — название групп, а номера — это абсолютные номера атомов, входящих в группу. htail — группа концевых углеродов лежащая «вверху» по оси Z, ltail — то же, но с другой стороны нанотрубки.
Синтаксис команды g_energy следующий:
$ g_energy -f *.edr -s *.edr -o ener.edr
При просьбе выбрать поля из таблице, необходимо выбрать полную потенциальную и кинтечискую энергии.
На выходе обоих программ получаются файлы XVG - текстовые файлы с таблицами интересующих значений, которые затем обрабатываются в среде IPython, позволяющей, в частности визуализировать изменения наблюдаемых параметров во времени (Рис.5).
Следующий этап работы: вычисление модуля Юнга на основе имеющихся таблиц энергий и координат. Вычисления производятся в MATLAB-подобной среде IPython. Для загрузки XVG-файла в массив используется функция load:
> a = load('file.xvg')
Для вычисления среднего значения массива — функция mean:
> eav = mean(a[:,3])
Для визуализации массива — функция plot:
13
> plot(a[:,0], a[:,1])
Вычисления модуля Юнга производятся согласно формуле 1. В одном случае в качестве силы F подставляется масса одного углеродного кольца, умноженного на ускорение, заданное в mdp-файле. Во втором случае в качестве F подставляется разница средних значений потенциальной энергии системы, деленная на среднее удлинение по оси Z.
Обратите внимание, что в обоих случаях необходимо брать среднее значение энергии и удлинения для равновесной области. Итерация, начиная с которой систему можно считать равновесной определяется качественно по графику энергии (см. Рис. 6).
Для вычислений необходимо учитывать размерности физических единиц, принятые в
GROMACS:
1. Энергия — кДж/моль
2. Ускорение — нм/пс
2 3. Расстояние — нм
4. Время — пс
Часть 3. Скорость диффузии воды в углеродной нанотрубке
Рассчитать коэффициент диффузии, исходя из достаточно коротких (200-500 пс) траекторий корректно невозможно из-за существенного вклада в суммарное смещение центров масс каждой конкретной молекулы случайных тепловых движений молекул.
Поэтому необходимо нормировать полученное для коэффициента диффузии воды внутри нанотрубки (Рис.7,8) значение, используя коэффициент диффузии воды в объемной фазе, так как этот параметр, при использовании более длинных (10-20 нс) траекторий хорошо согласуется с известным из эксперимента значением.
Для выполнения этой части работы необходимо сначала добавить в систему воду. Эту операция выполняется при помощи утилиты gen_box:
$ gen_box -sp *.gro -cs tip4p.gro -o *.gro
Как нетрудно видеть, для модели воды используется модель tip4p, в которую добавлена четвертая виртуальная частица, описывающая неподеленную электронную пару атома кислорода, что необходимо для корректного образования водородных связей в системе.
Файл молекулярно-динамических параметров подготавливается по уже известной схеме. Файлы TOP необходимо модифицировать в связи с добавлением в систему воды, а
14
SOL 256
.. в TOP-файл. После окончания молекулярно динамических расчетов необходимо извлечь из полученной траектории координаты молекул воды, которые необходимы для расчета коэффициента диффузии. Для этого используется утилита g_traj. Соответственно, ей необходимо передать номера атомов интересующих молекул воды, что делается при помощи соответствующего индекс-файла (*.ndx). Создание индекс-файла уже обсуждалось ранее. В данном случае необходимо внести группу молекул воды, находящихся внутри нанотрубки.
Затем из траектории извлекаются координаты атомов интересующих молекул воды:
$ g_traj -f *.trr -s *.tpr -n *.ndx -ox coord.xvg -mol
, где ключ [-mol] указывает на то, что будут выписаны не координаты атомов, а координаты центров масс молекул.
Обработка полученных данных производится также в среде IPython. Написать программу необходимо самостоятельно, руководствуясь формулой:
D = (X(t)-X(0))^2/2nt
.., где D – коэффициент диффузии, X(t) — радиус-вектор положения молекулы в момент t, n – размерность пространства. Параметр n будет отличаться для молекул воды в объемной фазе и внутри нанотрубки.
Требования к отчету
В конце практикума студенты готовят письменный отчет по задаче. Отчет должен включать следующие разделы:
−
Введение, включающее краткую характеристику метода молекулярной динамики, а также характеристика объекта исследования - нанотрубок, их структуры, свойств и возможных путей их практического использования.
−
Цель и основные задачи данного практикума.
−
Описание хода выполнения задачи, состоящее их трех частей, посвященных трем разделам задачи. Каждая часть должна начинаться с постановки задачи и краткого описания используемых для ее решения подходов. Затем должно следовать собственно описание хода работы, снабженное промежуточными числовыми
15
−
Заключение, в котором кратко формулируются полученные результаты и дается сравнительная характеристика нанотрубок различного строения.
Контрольные вопросы
1. Сколько атомов содержит минимальная кристаллографическая ячейка монослоя графита? Пояснить на схеме.
2. Каковы типичные размеры (радиус, длина) нанотрубки?
3. Рассчитать погрешности в определении равновесной длины связи в двухатомной молекуле и ее потенциальной энергии, при оптимизации с пороговым значением градиента энергии 10 кДж моль
-1
нм
-1
. Потенциальная энергия системы задана, как
U = k (R
AB
-R
0
)
2
/ 2.
4. Справедливы ли результаты гармонического анализа для системы не находящейся в потенциальном минимуме поверхности потенциальной энергии? Пояснить на примере одномерного потенциала.
5. Найти частоты и вектора нормальных колебаний одномерной системы состоящей из двух частиц массой m
1
и m
2
, если потенциальная энергия данной системы задана как
U = k
1
(x
1
-x
2
)
2
/ 2.
6. Что такое шаг интегрирования и радиус обрезания в методе молекулярной динамики?
16
Рис. 7: Цепочки связанных друг с другом водородными связями
молекул воды внутри нанотрубки.
Чем важны эти параметры?
7. Что такое периодические граничные условия в методе молекулярной динамики, для чего они используются?
8. В чем состоит метод управляемой молекулярной динамики?
9. Что такое модуль упругости? Какие существуют методы расчета коэффициента упругости на основании данных молекулярной динамики?
10. Что такое коэффициент диффузии? Какие микроскопические характеристики молекулярной системы его определяют?
11. Почему коэффициент диффузии воды внутри нанотрубки выше, чем коэффициент диффузии воды в объемной фазе?
Рекомендуемая литература
1. Сивухин Д.В. Общая физика. Москва. Наука. 1974.
2. Тихонов А.Н., Самарский А.А. Уравнения математической физики. Москва. Наука.
1972.
3. David van der Spoel, Lindahl E., Hess B., el. all. GROMACS User Manual. Department of
Biophysical Chemistry, University of Groningen. Nijenborgh 4, 9747 AG Groningen, The
17
Рис. 8: Молекулы воды внутри нанотрубки. Хорошо видно, что
молекулы воды располагаются внутри нанотрубки в весьма
ограниченном объеме.
Netherlands. 1991-2000.
4. Suzuki K. On Elastic Properties of Single-walled Carbon Nanotubes as Composite
Reinforcing Fillers. J. Comp. Mat. 41(9):1123-1135. 2007.
5. ChenXin F.U., YunFei C., JiWei J. Molecular dynamics simulation of the test of single-
walled carbon nanotubes under tensile loading. Sci China Ser E-Tech Sci. 50(1):1-17. 2007.
6. Garde A.S., Hummer G. Osmotic water transport through carbon nanotube membranes.
PNAS, 100(18):10175-10180. 2003.
7. Walthery J.H., Jaez R., Haliciogluz T., Koumoutsakosy P. Molecular dynamics simulations
of carbon nanotubes in water. Center for Turbulence Research Proceedings of the Summer
Program 2000 18