CC BY
65
Труды Карельского научного центра РАН № 1. 2013. С. 96-111
УДК 519.248.23:541.64:539.199
СУЩЕСТВЕННАЯ ВЫБОРКА ПРИ МОДЕЛИРОВАНИИ НЕПРЕРЫВНОГО СПЕКТРА КОНФОРМАЦИЙ МАКРОМОЛЕКУЛ МЕТОДОМ МОНТЕ-КАРЛО
А. Л. Рабинович1, Д. В. Журкин2
1 Институт биологии Карельского научного центра РАН,
2 Петрозаводский государственный университет
Описан алгоритм, разработанный для моделирования молекул цепного строения методом Монте-Карло в предположении о непрерывном спектре конформаций, с использованием процедуры существенной выборки. Алгоритм учитывает взаимозависимость каждых трех углов внутреннего вращения вдоль по цепи. Он использован для изучения равновесных свойств углеводородных цепей различного строения.
Ключевые слова: метод Монте-Карло, существенная выборка, цепные молекулы, непрерывный спектр конформаций.
A. L. Rabinovich, D. V. Zhurkin. AN IMPORTANCE SAMPLING FOR MONTE CARLO SIMULATION OF A CONTINUOUS SPECTRUM OF CONFORMATIONS OF MACROMOLECULES
An importance sampling technique for Monte Carlo simulation of chain molecules assuming a continuous spectrum of their conformations is described. The algorithm takes into account the interdependence of each three torsion angles along the chain. It is used to study the equilibrium properties of hydrocarbon chains of different structure.
Key words: Monte Carlo simulation, importance sampling, chain molecules, continuous spectrum of conformations.
Введение
Для изучения свойств молекулярных систем широко используются методы компьютерного моделирования [10]. Наиболее распространенными являются методы Монте-Карло (МК) и молекулярной динамики (МД). Первый из них позволяет получить средние значения макроскопических наблюдаемых величин по статистическому ансамблю, а второй -средние по времени. Однако используемые в методе МК выборки конечны, а доступные
в методе МД времена моделирования сравнительно невелики, поэтому важно проводить сравнение между собой данных, полученных в рамках различных подходов. Актуальной задачей является также совершенствование методов компьютерного моделирования. В настоящей работе описана специальная вычислительная схема метода МК, разработанная для генерирования конформаций цепных молекул, в предположении о непрерывном
изменении углов внутреннего вращения в полном диапазоне от 0 до 2п. Учтена взаимозависимость каждых трех углов внутреннего вращения вдоль по цепи. При генерировании использована существенная выборка.
Средние значения
Среднее значение < Н > любой макроскопической наблюдаемой величины Н(ж^ ж2, ..., жп) системы частиц (атомов данной молекулы) в каноническом ансамбле можно вычислить следующим образом [1, 11]:
< Н >=
/ / ... /п Н(жі, ж2, ..., жп) ■ ехр[ — и(жі, Ж2, ..., жп)/&£Т]^Ж1 ^Ж2...^ЖП
Здесь (Ж1,Ж2, ...,ЖП) странстве О этой системы, кв - постоянная Больцмана, Т - температура, и(ж1, ж2, ..., жп) -гамильтониан системы, в котором опущены члены кинетической энергии, поскольку декартовы координаты и импульсы частиц входят в гамильтониан (в потенциальную и кинетическую энергию) по отдельности; переменные разделяются, интеграл по всем импульсам берется тривиальным образом и его можно сократить в числителе и знаменателе. В итоге под фазовым пространством О подразумевается конфигурационная его часть, а под и(ж1, ж2, ..., жп) - потенциальная энергия, п -размерность пространства О. Если количество атомов в системе па, то п равно количеству декартовых координат (3па). Для молекул цепного строения, макромолекул конфигурационное пространство многомерно и функция и(ж1, ж2, ..., жп), как правило, очень сложна и в аналитическом виде не представлена. Поэтому провести интегрирование в выражении (1) и выразить < Н > в конечной аналитической форме в большинстве случаев не удается. Но можно использовать один из численных методов решения, - например, метод МК; литература о нем и его приложениях весьма обширна (см., напр., [1-3, 5-7, 10-12, 14, 15, 19, 25, 26]).
Метод Монте-Карло
В методе МК задачу сводят к моделированию случайных величин и статистической оценке их характеристик. В простейшей формулировке это означает следующее. Для того чтобы приближенно вычислить некоторую скалярную величину а, предлагается придумать такую случайную величину £, у которой существует математическое ожидание М£, равное а. Тогда для оценки величины а выбирается К независимых реализаций £1, £2,..., £к случайной величины £ и вычисляется среднее арифметическое £К = (1/К) • ^К=1 &. Последовательность одинаково распределенных независимых случайных величин, у которых существуют математические ожидания, под-
//.../п ежр[—и(ж^ ж2,..., жга)/йвT^ж^ж^^ж™ ( )
точка в фазовом про- чиняется закону больших чисел (см., напр.,
[4, 6, 12], теорема Хинчина), поэтому среднее арифметическое этих величин сходится по ве-
роятности к математическому ожиданию, т. е.
— р
при К — те имеет место £ к —— а. При больших К величина £к ~ а, и вместо а можно использовать [6, 12] оценку £К = (1/К) ^^=1 £*. По определению, последовательность случайных величин £1, £2,..., £к,... сходится по вероятности [12] к постоянной а, если для любого К > 0 вероятность Р{|£К — а| ^ К} — 0 при К — те (что равносильно Р{|£К — а| < К} — 1 при К — те). Здесь Р{|£К — а| ^ К} есть вероятность того, что £к удалена от а не менее чем на К. Смысл сходимости по вероятности в том, что вероятность значительного отклонения величины £к от величины а при большом К небольшая; тем не менее, сильно отклоняющиеся значения £к могут находиться сколь угодно далеко по К, но этих отклонений не очень много, вследствие чего вероятность их появления становится все меньше с ростом К.
Пусть ^£ - (конечная) дисперсия случайной величины £. Можно использовать понятие вероятной ошибки [12] гк = 0,6745 • (Я£/К )1/2, для которой Р{|£к — а| < гк} ^ Р{|£К — а| > Гк} ~ 1/2, т. е. одинаково вероятны ошибки, большие чем Гк, и ошибки, меньшие чем гк , и охарактеризовать порядок ошибки метода. Вероятная ошибка гк ~ К-1/2. Критерий вероятной ошибки гк полезен для сравнения степени "эффективности" [12] (или выгодности, удачности) разных алгоритмов, которые использованы (или можно использовать) для решения одной и той же задачи. Действительно, из выражения для Гк следует, что количества К осредняемых значений случайной величины К = ^£ • (0, 6745/гк)2. Пусть есть два алгоритма. Если задать вероятные ошибки гк для обоих одинаковыми, то К' = ££' • (0, 6745/гк)2, К'' = ££'' • (0, 6745/гк)2. Пусть £' и £'' - время, затрачиваемое на расчет одного значения случайных величин £' и £'', соответственно. Общее время расчета равно £' • К' = £' • ^£' • (0, 6745/гк)2
97
и Ь" • К'' = Ь" • ^£'' • (0, 6745/гк)2. Чем меньше время расчета (£' • К' или Ь'' • К''), которое необходимо для достижения заданной вероятной ошибки, при прочих равных условиях, тем более эффективным является алгоритм. А если ввести понятие "трудоемкости" [12] алгоритма, которое используется для сравнения алгоритмов и равно произведению дисперсии осредняемой случайной величины ^£ на время Ь расчета одного значения £ (т. е. Ь' • ^£' или Ь'' • ^£''), то можно сказать, что из двух алгоритмов более эффективным является менее трудоемкий.
Случайных величин £, - таких, что М£ = а, существует бесконечно много, и теория методов МК отвечает на вопросы о том [12], каким образом выбрать удобную величину £ для решения данной задачи и как находить значения £1, £2,... произвольной случайной величины £.
Вычисление интегралов
Представляется целесообразным кратко привести сведения о двух известных схемах применения метода МК для вычисления интегралов; они важны для последующего описания разработанного алгоритма.
Схема 1 ([12], с. 93). Пусть требуется приближенно вычислить интеграл ■] = /(ж,у) •
р(ж,у)^ж^у, где О - произвольная область плоскости ж, у, а р(ж, у) - некоторая плотность вероятностей: /ср(ж,у)^ж^у = 1. Для оценки .] методом МК вводится случайная точка Q в области О с плотностью вероятностей р(ж, у) и скалярная случайная величина ¥ = /(ж,у), -такая, что ее математическое ожидание М¥ = ■], т. е. равно искомому значению интеграла. Тогда, если Q1,...,Qк - независимые реализации случайной точки Q, а ¥1 = /^1),..., ¥к = / ^к), то согласно приведенной выше схеме оценкой .] будет величина [12] ¥ к =
(1/К) • ^К=1 ¥г. При этом ¥ К —— .
Важно отметить [12], что если О -ограниченная область, то по этой же схеме можно получить оценку и интеграла вида /с /(ж,у)^ж^у. Действительно, пусть площадь области О есть £с. Тогда плотность вероятностей случайной точки, равномерно распределенной в О, есть 1/5Ь. Если ввести /1(ж,у) =
• /(ж, у) и обозначить 1/£с = р1(ж,у), то этот интеграл приводится к предыдущему виду: /с /(ж,у)^у = /с /1 (ж, у) • Р1(ж,у)^ж^у.
Если искомую оценку получают усреднением по выборке значений функции с использованием плотности вероятностей случайной точки, распределенной в данной области равномерно (как в последнем примере), то такую
выборку в литературе по применению методов МК в физике (статистической, теоретической) принято называть ’’simple” sampling [15, 25], т. е. "простая" выборка [2] (иногда используется термин "прямая" выборка [14]).
Схема 2 ([12], с. 108). Пусть J0 =
Sg f (x, y)dxdy, где G - ограниченная или неограниченная область, и J'g |f(x,y)|dxdy > 0. В области G вводится плотность вероятности р(ж,у). Если р(ж,у) > 0 в тех точках, в которых f(ж, y) = 0, то р(ж, y) называется допустимой по отношению к f (ж, y); допустимая плотность р(ж, y) может обращаться в нуль только в тех точках, где f (ж, y) = 0. Если же р(ж, y) > 0 всюду в G, то она допустима по отношению к любым f (x,y). Пусть Go - множество точек, в которых f (ж, y) = 0, и пусть G+ = G — G0. Вводится функция F0(x, y), следующим образом: F0(x,y) = f (ж, y)/p(x, y), если точка (ж,у) принадлежит области G+, и F0(ж,у) = 0, если точка (ж, у) принадлежит области G0. Здесь р(ж,у) - допустимая плотность. Для оценки J0 методом МК, как и в схеме 1, вводится Q - случайная точка, определенная в области G с плотностью вероятностей р(ж,у). Математическое ожидание величины F0(Q) есть MF0(Q) = J'G F0(ж,у) ■ р(ж,у)^ж^у = /g+ f (ж,у)ЛЫу = J0 ввиду того, что функция f(ж,у) =0 в точках (ж,у), не принадлежащих области G+. Итак, MF0(Q) равно искомому значению интеграла. Дальнейшие действия соответствуют схеме 1: если Qi,...,Qk -независимые реализации случайной точки Q,
а Fi0 = f(Qi)MQi),..., FK = fQk)/p(Qk),
то оценкой J0 будет величина F0k = (1/K) ■
F0
Z^i=1 .
Очевидно, плотность р(ж, y) желательно выбрать так, чтобы минимизировать дисперсию DF0(Q). Доказано [12], что минимальная дисперсия DF0(Q) реализуется в случае, когда плотность р(ж,у) пропорциональна |f(ж, у) |; при этом она равна |f(ж, у) | ■ [/G |f(ж,у)|йж^у]-1. Ясно, однако, что несмотря на привлекательность достижения минимальной дисперсии, строго следовать этой ре-комендгции нельзя, поскольку задача вычисления Jg |f (ж,у)|^ж^у по трудности эквивалентна исходной задаче вычисления искомого интеграла. Но можно сделать и другой вывод: для получения величин F0 с небольшими дисперсиями желательно выбирать плотность р(ж, у) по возможности пропорциональную |f (ж, у) | [12], - тогда в тех частях области G, в которых |f(ж, у)| больше и вклад которых в J0 более существен, будет выбираться больше случайных точек [3, 12]. Иными словами,
чем меньше меняется функция f (ж,у)/р(ж,у) (чем ближе она к постоянной), тем меньше окажется дисперсия DF0(Q). В этом состоит основная ценность данного приема; он предложен в работах [22, 23] и широко применяется [1, 2, 7, 12, 14]. Такая выборка называется "importance" sampling [15, 25] или "существенная" выборка [3, 7, 12], выборка "по значимости" [2, 5], выборка "по важности" [1, 6, 7], "предпочтительная" выборка [14].
Общий подход к оценке средних
Обе схемы, описанные выше для двух переменных, обобщаются на случай многих переменных. Возвращаясь к исходному выражению (1) для среднего значения < H >, перейдем от точек (ж, у) плоскости к n-мерным точкам (ж1, ж2,..., жп) конфигурационного пространства Q. Интеграл в знаменателе выражения (1) есть так называемый
статистический интеграл 2 (для дискретного случая - статистическая сумма [9]), 2 = ... п ежр[—и(ж1,..., жга)/йвТ]^ж1...^жга, а величина (1/2) ■ ежр[—и(ж1,..., жга)/йвТ] есть плотность вероятности точки (ж1 , ...,жп). Для приближенного вычисления интегралов в выражении (1) можно применить как "простую", так и "существенную" выборку.
Реализуем "простую" выборку. п-мерная точка (ж1,...,жп) является аналогом рассмотренной выше 2-мерной точки ^. Независимыми реализациями ^ случайной точки Q будут точки (ж1,..., жП), где і = 1, 2,..., К. Пусть они распределены равномерно в О. Скалярной случайной величиной ¥ (из схемы 1) является выражение (1/2) ■ Н(ж1,...,жп) ■ ежр[—и(ж1, ...,жга)/&£Т], и его математическое ожидание, согласно (1), равно < Н >. Оценка Нк выражения (1) методом МК будет иметь вид [1, 2, 11]
-П Ei=i H(xi,x*2,..., xn) ■ exp[—U(xi,x*2,..., xn)^T]
H к = --------------T?------------:---:-----------------------
Ei=i exp[—U(xl, x2,..., xn)^T]
— р
При этом Нк —— < Н >. Отметим, что выражение (2), в силу специфики исходного выражения (1), имеет вид формулы усреднения в каноническом ансамбле для случая, когда переменные имеют только дискретный спектр. Однако в данном случае суммирование производится не по дискретному спектру (всем состояниям), а по случайным точкам непрерывного спектра, и К является не количеством состояний, а размером выборки. В выражении
(2) множитель ежр[—и(ж^, ..., ж^)/квТ] может изменяться на многие порядки и, следовательно, более предпочтительной была бы "существенная" выборка.
Вариант "существенной" выборки, который можно было бы назвать наиболее эффективным, выгодным, т. е. наименее трудоемким (в том смысле, который описан выше), очевидно, был бы достигнут в случае использования функции (1/2)• ежр[—и(ж1, ж2, ..., Жп)/квТ] в качестве плотности распределения случайной точки в области О. Пусть точки (ж1, ж2, ..., жП) выбираются в О не равномерно, как в выражении (2), а именно с такой плотностью. Тогда оценка выражения
(1) сводится к среднему арифметическому
НК = (1/К) • ЕК=1 Н(ж1,ж2,...,жП). Несмотря на столь простой вид итоговой оценки, реализовать эту идею невозможно, поскольку по трудности задача предварительного вычисления этой плотности сопоставима с задачей вычисления исходного выражения (1): речь идет о необходимости предварительного вычисления потенциальной энергии системы частиц и(ж1, ж2, ..., жп). Тогда нужно подобрать такую плотность р(ж1, ж2, ..., жп), которая была бы допустимой по отношению к подынтегральному выражению в (1) и, насколько возможно, пропорциональной ему, но при этом могла быть рассчитана предварительно, до проведения расчета методом МК. В итоге можно добиться дисперсии, значительно меньшей, чем в алгоритме простой выборки, и более эффективного (менее трудоемкого) алгоритма. Пусть некоторую плотность р(ж1, ж2, ..., жп), удовлетворяющую этим условиям, найти удалось. Повторяя действия, описанные в схеме 2 выше, получим следующую оценку [1, 2, 11, 14]:
Т7 Ei=i H(xl, x2,..., xn) ■ exp[—U(xi, x2,..., xn)/kвT]/p(xi, x2,..., xn)
H к = ------------к-----------•--•-----------------•--•---------------
Ei=i exp[—U(xi, x2,..., xn)/kвT]/p(xi, x2,..., xn) --------------------------------------------------------(99)
Различие между (3) и (2) обусловлено именно способом выбора точек в конфигурационном пространстве.
Если выбрать р(ж1, ж2, ..., жп) в (3) пропорционально ежр[—и(ж1,ж2, ...,Жп)/квТ], т. е. заранее вычислить какую-то часть энергии и(Ж1, Ж2, ..., жп), то получим искомый результат. Действительно, энергию и можно представить, например, в виде суммы двух слагаемых (назовем их энергией ближних и^ог* и энергией дальних и^ога^ взаимодействий), и = и^ог* + и^^. Тогда ежр[—и(ж1, ж2, ..., Жп)/квТ] будет произведением ежр[—и^ог*(Ж1, Ж2, ..., Жга)/квТ]
и ежр[—и^огад(Ж1, Ж2,..., жп)/квТ]. Выделение и5^ог* из общей энергии и следует осуществить таким образом, чтобы была возможность вычислить ее на предварительном этапе,
до осуществления процедуры генерирования конформаций молекул методом МК. Детали реализации этой процедуры будут описаны ниже. Если и^ог* вычислена, то в качестве плотности р(ж1 , ж2,...,жп) распределения случайных точек можно использовать (1/2*ог*) ■
вЖ^р[ ийЬог*(ж 1, Ж2, ..., Жп)/kBT], где 2«Ь,ог* —
/ / ... /п ежр[ — и^ог*(Ж1, ...,Жп)/квТ]^Ж1...^Ж„. Ясно, что это выражение является допустимой плотностью по отношению к подынтегральной функции: оно ^ 0 во всей области конфигурационного пространства, а если обращается в нуль при и^ог* — те, то в этих же точках обращается в нуль и экспонента, в которую входит полная потенциальная энергия системы и, поскольку тогда и — те. В итоге вместо
(3) получим следующую оценку
— _ Ег=1 Н(ж*1 ,ж2,...,жП) ■ ежр[—и^опд(ж*1 ,ж2, ...,жП)/квТ]
К К ’ ’ ■
Ег=1 ежр [ — и^опд (ж1 , Ж2, ..., ЖП ) /квТ]
При компьютерной реализации эту плотность удается воспроизвести лишь приблизительно, и итоговые формулы (см. ниже) отличаются от (4). Тем не менее, выборка оказывается в итоге "существенной" .
Отметим, что при изучении свойств макромолекул идея вычисления или аппроксимации части энергии и(ж1, ж2, ..., жп) в методе МК использовалась многократно, предложено много вариантов для разных молекулярных систем, и в каждом из них учтена специфика этих систем и характер задачи; ряд вариантов упомянут, например, в обзоре [11]. В настоящей работе эта идея использована для разработки алгоритма, который можно применить для исследования свойств молекул цепного строения определенных классов, и в нем учтена специфика этих цепей. В частности, ниже алгоритм расчета реализован для углеводородных цепей вида СНз — (СН2)/ — (СН = СН — —
(СН2)Й — СНз, где /, ^,д - целые.
Описание цепных молекул
Приведем основные сведения о принятом описании цепных молекул. Конформации мак-
ромолекулы могут быть охарактеризованы либо декартовыми координатами всех атомов, либо внутренними переменными: совокупностями валентных связей {1}, валентных углов {0}, углов внутреннего вращения {у} (т. е. торсионных углов). В рамках классической статистической механики существуют две схемы приближений [20], - "классическая гибкая модель" (согласно которой макромолекула - это система материальных точек, взаимодействие которых описывается некоторыми потенциалами) и "классическая жесткая модель" (согласно которой потенциалы действуют лишь между определенными участками макромолекулы, остальные же участки скреплены между собой геометрическими связями).
"Гибкая" модель при изучении равновесных свойств считается более предпочтительной [18, 20] (несмотря на ряд нерешенных вопросов [18, 20]), и она позволяет для среднего значения < Н > вместо выражения (1) в декартовых координатах записать выражение, в котором все величины зависят от внутренних переменных:
Н > Я... / Н ({у}, {0}, Щ) ■ ежр[—и ({у}, {0}, М)/кв ТИуДОИЩ
Я ..., ежр[—и ({у}, {0}, Щ)/квТИу}<({0}<ВД ■ ( )
Здесь, как и в (1), Н - некоторая макроскопи- внутреннего вращения, валентных углов и ва-
ческая наблюдаемая величина изучаемой си- лентных связей, соответственно, которые об-
стемы частиц (молекулы), и - потенциальная разуют конфигурационное пространство.
энергия; {у}, {0} и {1} - совокупности углов
Для получения оценки величины (5), аналогичной (3), методом МК можно использовать существенную выборку. Пусть случайные точки в конфигурационном пространстве есть {у}^, {0}^ и {1}^, а количество конформаций (состояний V) молекулы в выборке есть
ш (для внутренних переменных обозначение
V использовано вместо г, а ш - вместо объема выборки К). Тогда, выбирая случайные точки с некоторой плотностью вероятностей р({у}, {0}, {1}), получим следующую оценку Н ш:
Н _ Е"=1 н({уГ, {0Г, ОГ) ■ ехр[-и({уГ, {0Г, ОГ)Авт]/р({уГ, {0^, {^) (б)
" Е"=1 ехр[-и({у}-, {0}-, {1}-)/кБТ]/р({у}-, {0}-, {1}-) • ( )
Выражение (6) можно сделать более компакт- менных. Обозначая Н({у}^, {0}^, {1}^) _ Ни; ным, если убрать явное перечисление совокуп- и({у}^, {0}^, {1}^) _ ии; р({у}^, {0}^, {1}^) _ ностей случайных значений внутренних пере- , приведем (6) к виду
Нш _
Силовое поле
Рассмотрим вопрос об общем подходе к расчету энергии цепных молекул. Как правило, в расчетах используется метод атом-атомных потенциальных функций. В литературе существуют наборы таких функций и апробированные параметры (силовые поля) для цепных молекул разных типов. В частности, для молекул липидов и их компонентов - углеводородных цепей рассматриваемого типа предложено поле СНЛИММ27 [16, 17, 24, 28]; в настоящей работе использован модифицированный вариант этого поля [21]. Энергия и системы частиц вычисляется в виде суммы энергий ва-
N
лентных связей и _ Е Кгп ■ (1„ - 1о„ )2 и уг-
П=1
N.
лов и _ ^ Кв„ ■ (0п - 0о„)2, энергии Юри-
П=1
N
Брэдли ии _ ^ Кип ■ («п - ио„)2, энергии
п=1
Nn
неплоских отклонений атомов ич _ Е КЧп ■
п=1
N
(Пп - Поп)2, торсионной энергии и^ _ ^ К^п ■
п=1
[1 + С08(рга ■ уп - £га)], энергии невалентных взаимодействий игаЬ _ Е в, ■ [(Дт»п,3/Гу)12 -
г,,
2 ■ (Ятт^/г,)6] и электростатической энергии
ие1 _ ЕК^ ■ 9,-)/(ве1 ■ Гг,)]. г,,
Здесь N1, N0, Жи, N, - соответственно,
количество валентных связей; валентных углов; различных пар атомов, разделенных двумя валентными связями; углов неплоских отклонений в молекуле; торсионных углов. К^, К#п, Кип, КЧп, К^п - силовые постоянные для
(7)
расчета перечисленных компонентов энергии, соответственно. Среди параметров слагаемых энергии используются текущие и равновесные значения: длины п-й валентной связи (1п и 10п), п-го валентного угла (0п и 0Оп), п-го расстояния между атомами, разделенными двумя валентными связями (ип и иоп), п-го угла неплоских отклонений (пп и Поп). уп - это значение п-го торсионного угла, а рп, - мно-
житель и сдвиг фаз; г, - расстояние между атомами с номерами г и ], разделенными более чем двумя валентными связями; в, и Ятт„ -глубина и положение соответствующего минимума энергии взаимодействия; дг, 9, - парциальные заряды на атомах, - диэлектрическая постоянная.
Сравнение внутренних переменных
Переменные {у}, {0}, {1} не являются, в определенном смысле, равноправными между собой. Действительно, при данной температуре, в условиях, близких к состоянию равновесия молекулы, каждый угол внутреннего вращения из совокупности {у} может с высокой вероятностью изменяться во всем диапазоне [0, 2п], тогда как изменения валентных связей {1} и валентных углов {0} происходят лишь в малых областях вокруг их равновесных значений (это так называемые "малые колебания"). Данное различие использовано при разработке настоящей вычислительной схемы метода МК: генерирование каждой конформации макромолекулы осуществляли в 2 этапа.
На 1-м этапе генерирования конформации
V все величины {1} и величины {0} (или их большинство) фиксировали при равновесных значениях, т. е. соответствующие компонен-
Е"=1 Ни ■ ехр[-и*/квТ]/р^ Е"=1 ехр[-и -/кв Т]/р-
■©
ты энергии обращали в нуль. Конформацию номер V генерировали только по торсионным углам {у}. Последние изменялись в диапазонах [0,2п], но если рассматриваемая молекула содержала двойные связи е1э (вращения вокруг которых сопряжены с преодолением большого потенциального барьера), то, аналогично валентным углам, соответствующие углы вращения у фиксировали в положении равновесия. При генерировании конформации V по совокупности углов {у} был предложен вариант существенной выборки, детали которого будут описаны ниже. В отличие от углов внутреннего вращения, независимое варьирование в методе МК всех валентных углов осуществить невозможно, - например, углов, имеющих вершину на одном и том же атоме. Поэтому на 2-м этапе осуществляли переход в полученной конформации от внутренних переменных (валентных связей, валентных углов, торсионных углов) к декартовым координатам атомов и проводили варьирование всех координат, - независимо для каждого атома. При этом было учтено, что, поскольку валентные связи {1} и валентные углы {0} претерпевают лишь малые колебания вокруг положений равновесия, изменения декартовых координат атомов также
могут быть лишь малыми. Оценка величины диапазона варьирования декартовых координат проведена из условия, чтобы максимальные изменения значений энергии валентных углов и# и валентных связей и^, отвечающих этому диапазону, в сумме не превышали изменений торсионной энергии цепи и энергии невалентных взаимодействий ипь. Расчеты показали, что этому условию отвечает диапазон ± 0,02 А вокруг равновесных положений всех атомов. Ввиду малости этого диапазона полагали, что в его пределах значения координат распределены равномерно; процедуру генерирования координат осуществляли с использованием генератора (датчика) псевдослучайных чисел [27], равномерно распределенных в интервале [0,1).
Итак, валентные связи I и валентные углы в сначала фиксировали при равновесных значениях. В выражении (5) можно явно указать поэтому только торсионные углы. Пусть N - количество атомов основной цепи молекулы, У1,У2...,УМ-1 - углы внутреннего вращения вокруг связей основной цепи, и(у1, у2,..., у^-1) - потенциальная (конфор-мационная) энергия. Тогда вместо (5) будем иметь
< Н >=
2п 2п 2п
/ / ... / Н(У1,У2, ...,ум-1) ■ ехр[—и(У1,У2, ...,ум-1)/квТ]^У1^У2...^Ум-1 0 0 0
2п 2п 2п
/ / ... / ехр[ —и(У1,У2, ...,Ум-1 )/квТ]^У1 ^У2...^Ум-1
0 0 0
Энергия и(у1, У2,..., Ум-1) при рассмотрении системы макромолекула - растворитель является энергией системы, состоящей из макромолекулы и окружающих молекул растворителя, усредненной (при любом данном состоянии У1,У2...,Ум-1 макромолекулы) по всем поло-
жениям и ориентациям молекул растворителя [18].
При использовании существенной выборки, выбирая случайные точки с некоторой допустимой плотностью р(у1, У2, ..., Ум-1), получим следующую оценку Нш методом МК:
~тт ^=1
Н ш
Е Н(У1, У2,..., Ум-1) ■ ехР[-и(У1, У2,..., Ум-1 )/квТ]/Р(У1, У2,..., Ум-1)
Е ехр[-и(у 1, у2,..., ум-1)/квТ]/р(у1 , у2,..., уМ-1) ^=1
(9)
Расчет энергии цепной молекулы
Опишем более подробно процедуру расчета энергии и, а затем выбора допустимой плотности р(у1, у2,..., ум-1). Как упоминалось в разделе «Общий подход к оценке средних», потенциальная энергия всей молекулы и была
представлена суммой энергий ближних и*ог* и дальних и^ога^ взаимодействий, и = и^ог* + и^ога^. Выделение и^ог* проведено так, чтобы ее можно было вычислить на предварительном этапе, до осуществления процедуры генерирования конформаций методом МК.
Изучаемая цепная молекула была условно разделена на небольшие фрагменты, - такие, что каждый из них содержал 3 последовательных угла внутреннего вращения вокруг связей основной цепи, у7, у7+1, у7+2, варьирование которых приводило к изменению взаимных положений атомов. При конструировании цепи принято следующее правило: фрагменты соединяются не "встык", а с наложением друг на друга, - так, чтобы последовательные (соседние, смежные) фрагменты имели по два общих варьируемых торсионных угла из трех. Этим учтена непрерывная взаимозависимость каждых трех углов вдоль по цепи.
Энергия и8Ьог* является суммой энергий иШт этих фрагментов цепи, а энергия взаимодействий более удаленных вдоль по цепи атомов входит в игоп^. В итоге и8Ьог* = м-3
]Т иШт(у7,у7+1 ,у7+2). Здесь т7 - индекс
7=1
фрагмента: т - его идентификационный номер, который характеризует химическое строение фрагмента (т = 1, 2,..., т/, где т/ - общее количество разных фрагментов, необходимых для конструирования избранных молекул по указанному правилу), а 7 - его номер в последовательности фрагментов вдоль по цепи. Номера т для некоторых 7 в индексах т7 в сумме и8Ьог* могут быть одинаковыми. Если молекула - это полимерная цепь, то номера т будут циклически повторяться для одинаковых мономерных единиц.
Как упоминалось, алгоритм расчета реализован для углеводородных цепей вида СН3 — (СН2) / — (СН = СН — СН2 )* — (СН2 )* — СНз, где /, й, д - целые. Анализ показал, что любая из цепей такого вида может быть сконструирована должным выбором нескольких фрагментов из следующей совокупности (три связи, углы вращения вокруг которых варьируются, ниже выделены дефисом):
СН3 — СН2 — СН2 — СН2СН2;
СН3СН2 — СН2 — СН2 — СН2СН2;
СН2СН2 — СН2 — СН2 — СН2СН2;
СН2СН2 — СН2 — СН2 — СН2СН;
СН2СН2 — СН2 — СН2 — СНСНСН2;
СН2СН2 — СН2 — СНСН — СН2 СН2; СН3СН2 — СН2 — СН2 — СН2СН;
СН2СН2 — СН2 — СНСН — СН2 СН; СН2СН2 — СНСН — СН2 — СНСНСН2;
СН3 — СН2 — СНСН — СН2СН;
СН3СН2 — СНСН — СН2 — СНСНСН2; СН2СНСН — СН2 — СНСН — СН2 СН; СН2СНСН — СН2 — СН2 — СН2 СН3; СНСН2 — СН2 — СН2 — СН3;
СНСН2 — СНСН — СН2 — СН2 СН3; СН2СНСН — СН2 — СН2 — СН3.
Один из этих фрагментов представлен в виде схемы на рис. 1.
н н н н н н
н н н н н н
Рис. 1. Охема фрагмента СН2СН2 — СН2 — СНСН — СН2 СН2. Три варьируемых торсионных угла (у1, у2, уз) обозначены стрелками
Весовые множители
Энергию Um (^1, ^2, ) каждого из 16
представленных выше фрагментов вычисляли (с учетом только компонентов U<^, Unb и Uei) таким образом, чтобы сумма энергий любой совокупности фрагментов оказалась равной энергии ближних взаимодействий Ushort той цепи, которая составлена из этих фрагментов. Для этого при вычислении каждой Um вводили для трех перечисленных компонентов энергии весовые множители (1/3, 1/2 или 1). Множители были согласованы между смежными фрагментами и исключали возможность того, что одни и те же слагаемые окажутся дважды (или более раз) просуммированы в итоговой энергии Ushort всей цепи. Избран принцип равноправного относительно смежных фрагментов учета каждого компонента энергии. А именно, слагаемые энергии, зависящие строго от одного торсионного угла, учитывали с множителем 1/3 в трех фрагментах (вдоль по цепи), в которых этот угол варьируется. В первом фрагменте из этих трех угол варьируется как ^3, во втором смежном фрагменте - как ^2, в третьем смежном фрагменте - как ^1. Слагаемые, зависящие строго от двух торсионных углов, учитывали с множителем 1/2 в каждом из двух смежных фрагментов, в которых эти углы варьируются (в первом фрагменте - как ^2, ^э, а во втором смежном фрагменте - как ^1, ^2). Слагаемые, зависящие строго от трех торсионных углов, учитывали с множителем 1 только в том фрагменте, в котором эти углы варьируются. В энергии Um каждого фрагмента из рассмотрения были исключены такие слагаемые, которые формально зависят от трех (или двух) углов вращения, но в данном фрагменте m один из них не варьируется. Эти слагаемые включали в энергию Um того смежного фрагмента,
-(ш)
в котором варьируются все три (или, соответственно, два) торсионных угла.
В концевых фрагментах весовые множители некоторых компонентов энергии отличались от таковых для "срединных" фрагментов (для корректного расчета и*ог* всей цепи). Например, в концевых фрагментах все слагаемые энергии, которые зависели только от одного торсионного угла, учитывали полностью (с множителем 1, а не 1/3), если этот угол был первым (у1), поскольку ни в один из последующих фрагментов эти слагаемые больше не входят. Все слагаемые энергии, которые в концевых фрагментах зависели от второго угла (у2), учитывали не с множителем 1/3, ас 1/2, поскольку эти слагаемые входят не в три, а в два фрагмента: в концевой и смежный с ним фрагмент. Множители 1/3 для слагаемых энергии в концевых фрагментах были использованы только для третьего угла (у3). Табулирование значений энергии ит(у1 ,у2,у3) по трем углам вращения осуществлено с интервалом 1°.
Существенная выборка
Если энергия и^ог* исследуемой цепной молекулы рассчитана (по фрагментам), то алгоритм существенной выборки при генерировании конформаций этой молекулы методом МК будет эффективным, если организовать выборку совокупности случайных значений торсионных углов с плотностью вероятностей р — (1/^«^ог*) • ежр[ и«^ог*(у1 ,у2,...,уМ— 1)/квТ].
2п 2п
Здесь ^зЬог£ — / ... ]‘ ежр[ Ushoгt(у1, ...,
0 0
ум—1)/квТ] • ^уь..^ум —1.
М—3
Поскольку и^^ог* = ит7 (у7,у7+1,у7+2),
7=1
то числитель
ежр[ и^^ог*(у1, ..., ум—1)/квТ]/
м—3
— П ежр[—итт(у7,у7+1,у7+2)/квТ],
7=1
а
2п 2пм—3
/ ... I П ежр[ ит7 (у7, у7+1 ,у7+2)/ 0 0 7=1
кв Т] • ^у1...^ум—1.
Методика реализации этой идеи такова. Используя табулированные значения энергии ит каждого фрагмента, при данной температуре
Т рассчитывали, с интервалом 1°, функции ежр[—ит(у1 ,у2,у3)/квТ] и численными методами вычисляли интегралы
2п 2п 2п
/ / / ежр[—ит(у1,у2,у3)/квТ]^у1^у2^у3.
000
Пусть величина такого интеграла для фрагмента т равна /т. Для каждого т области от 0 до 2п по каждому из трех углов (у1,у2,у3), т. е. "куб" , разбивали на 1003 — 1000000 параллелепипедов (со-
стояний) так, чтобы вероятности осуществления каждого из состояний были равны друг другу. Указанное разбиение осуществлено последовательно, численно, - с использованием метода прямоугольников [8]. Сначала были рассчитаны такие значения углов по оси у1 (положения плоскостей), которые разбивают трехмерное пространство (у1,у2,у3) на 100 "слоев" , интегралы по которым от функции ежр[—ит(у1 ,у2,у3)/квТ] равны /т/100. Каждый из слоев затем разделили по углу у2 на 100 равновероятных "подслоев" из условия, чтобы интегралы от функции ежр[—ит(у1, у2, у3)/квТ] по каждому подслою были равны /т/10000, т. е. рассчитали значения углов по оси у2, которые разделяют между собой подслои в каждом слое. Наконец, каждый подслой разделили по углу у3 на 100 состояний (параллелепипедов) так, чтобы интегралы от функции ежр[—ит(у1 ,у2,у3)/квТ] по каждому параллелепипеду были равны /т/1000000 (рассчитали соответствующие граничные значения углов между параллелепипедами по оси у3 в каждом подслое). Очевидно, что при таком способе разбиения границы параллелепипедов в областях минимумов энергии каждого фрагмента т сгущаются, т. е. количество состояний -параллелепипедов в областях вблизи минимумов энергии оказывается больше, чем в областях вблизи максимумов. Массив чисел -границ между 1000000 параллелепипедов по всем трем углам у1, у2, у3 данного фрагмента используется далее для генерирования состояний с плотностью вероятностей, пропорциональной величине ежр[—ит(у1, у2, у3)/квТ], которая входит в качестве сомножителя в плотность р.
Например, на рис. 2 схематично приведено несколько плоскостей, которые делят область изменения трех варьируемых углов вращения ("куб" ) некоторого фрагмента на слои по оси у1, - это слои вблизи середины (180°) диапазона изменения угла у1 .
©'
Рис. 2. Разбиение части области изменения трех варьируемых углов любого фрагмента на равновероятные слои по оси у1
Зависимость энергии фрагмента, изображенного на рис. 1, от углов у2 и у3 при у1 — 180°, т. е. в одном из таких слоев, представлена в виде поверхности на рис. 3:
Рис. 3. Зависимость энергии (ккал/моль) фрагмента рис. 1 от углов ^2 и при = 180°. Отсчет энергии - от глобального минимума энергии данного фрагмента
Эту поверхность можно представить в виде эквиэнергетических линий, т. е. конформаци-онной карты, рис. 4:
Рис. 4. Конформационная карта фрагмента рис. 1 при у1 = 180°. Числа у кривых - энергия в ккал/моль, отсчет энергии - от глобального минимума
Рис. 5. Проекция границ параллелепипедов на плоскость (у2, у3) при у1 = 180° фрагмента рис. 1, совмещенная с соответствующей конформацион-ной картой. Числа у кривых - энергия, ккал/моль
Слой, в который попадает значение у1 = 180°, как и другие слои, разделен на подслои по оси у2 (’столбики’), а каждый подслой -на параллелепипеды по оси у3. Проекция границ итоговых параллелепипедов на плоскость у2, у3, в которой изображена также конфор-мационная карта, имеет вид, представленный на рис. 5.
Проекция границ параллелепипедов плоскость ^2, представлена на рис. 6:
на
105
О 45 90 135 180 225 270 315 360
Ч>2
Рис. 6. Проекция границ параллелепипедов на плоскость (у2, Уз) при у1 = 180° для фрагмента, изображенного на рис. 1
Легко видеть, что совокупность проекций границ параллелепипедов на рис. 6 воспроизводит очертания конформационной карты. Эти границы делят ось у2 фрагмента на полосы неравной ширины (это проекции "подслоев" ), а каждую из полос по оси уз — на прямоугольники (это проекции параллелепипедов) неравного размера. По количеству прямоугольников на единицу площади можно определить положения областей, в которых расположены минимумы энергии, - в них границы прямоугольников, действительно, сгущаются.
Необходимо также отметить следующее. Для расчета свойств цепных молекул часто применяется поворотно-изомерная модель [13], в которой используются только дискретные состояния, отвечающие положениям минимумов энергии. Из рис. 3, 4 видно, что потенциальная поверхность данного молекулярного фрагмента сложна и имеет пологие, асимметричные минимумы. Такая ситуация характерна и для потенциальных поверхностей при других значениях угла у1 данного фрагмента, а также для потенциальных поверхностей других фрагментов молекул данного типа. Если ярко выраженных, симметричных минимумов энергии нет, то использование поворотно-изомерного приближения затруднено, поскольку открытым остается весьма существенный вопрос о выборе значений параметров [13]. При использовании в настоящей работе непрерывного спектра конформаций учитывается асимметрия заселенности состояний. Это позволяет более адекватно моделировать макромолекулы заданной структуры, обеспечивает более широкие возмож-
ности в описании их характеристик в разных условиях, а также предсказании свойств в тех случаях, когда экспериментальные данные отсутствуют. Это создает основу для более корректного рассмотрения, например, эффектов, в которых важную роль играет "накопление" вдоль по цепи небольших отклонений углов внутреннего вращения от поворотноизомерных состояний. В этом состоят преимущества использования непрерывного спектра конформаций цепных молекул в теоретических моделях и, в частности, целесообразность разработки настоящего алгоритма. В последнем учитываются и другие важные факторы: воспроизводится точное химическое строение цепи и учитывается взаимозависимость каждых трех последовательных углов внутреннего вращения.
Генерирование конформаций
При генерировании конформаций цепи с углами у1,у2,уз, ...,ум-1 составляли надлежащую последовательность фрагментов, воспроизводящих строение данной цепи, согласно принятому правилу: любые два соседних (смежных) фрагмента должны иметь два общих варьируемых торсионных угла из трех. Номерам каждой тройки углов во фрагментах присваивали последовательные номера
(у1,у2,уз), (У2,Уз,У4), (Уз,У4,Уб) и т. д. вдоль по цепи. Для выбора значений первой тройки углов (у1,у2,уз) в данной цепи методом МК выбирали один из 1000000 параллелепипедов из надлежащего концевого фрагмента, т. е. номера трех сторон параллелепипеда: номер к(1) слоя (от 1 до 100), номер к(2) подслоя (от 1 до 100) и номер к(з) состояния в подслое (от 1 до 100). Каждый номер выбирался с помощью генератора псевдослучайных чисел [27], равномерно распределенных в интервале [0,1). Пусть € [0,1) - такие числа (г=1, 2, 3), тогда к® = ШТ(100 ■ пМ) + 1, где ШТ -целая часть числа.
В пределах избранного параллелепипеда выбирали фиксированные значения углов У1, у2 и уз, полагая, что ввиду большого количества (1000000) параллелепипедов, на которые разделено конфигурационное пространство каждых трех углов, изменениями плотности вероятности внутри каждого параллелепипеда можно пренебречь, т. е. распределение значений углов в нем считали равномерным (отметим, что при этом также вычисляли и реализуемую плотность вероятностей генерирования углов; см. ниже). При выборе этих значений использовали информацию о положе-
нии границ по осям y1, y2 и уз уже избранного параллелепипеда и три очередных псевдослучайных числа n(i) £ [0,1); i=4, 5, 6. Например, если y1fc(i), yi(fc(i)+i} - значения границ k(1) -го слоя по углу yi, то фиксированное значение
(fix)
y1k(1) угла yi внутри этого слоя вычисляли согласно y(1'kix))=уШ1)+(y1(fc(i)+1)- y1fc(D) ■ п(4). Аналогично, для расчета фиксированного значения f угла y2 внутри k(2) -го слоя использовали положения соответствующих гра-
(5) (fix)
ниц и число п , а для y3fc(3) угла y3 внутри
к(3)-го слоя - положения границ и число п(6).
Очередной угол y4 выбирали из разбиения, отвечающего смежному фрагменту с тройкой углов y2,y3,y4. Но при этом рассматривали лишь те 100 параллелепипедов смежного фрагмента, в которые попали уже избранные на предыдущем шаге фиксированные значения двух предыдущих углов, у>2 и у>3. Аналогично вышеприведенной процедуре выбирали некоторый номер k(4) одного из 100 этих параллелепипедов, для чего генерировали очередное (7-е) псевдослучайное число п(7) £
[0,1) и вычисляли k(4) = INT(100 ■ п(7) + 1). Затем внутри избранного параллелепипеда выбирали фиксированное значение угла y4. По-
следняя процедура, как и в приведенном выше примере, выполнялась с использованием известных значений границ по оси у4 данного параллелепипеда и еще одного генерированного псевдослучайного числа € [0,1).
Процедуру выбора фиксированных значений всех последующих углов у7, где 7 = 5, 6, ..., N — 1, проводили аналогично процедуре выбора значения угла у4 (выбирали номер одного параллелепипеда из 100 и затем - фиксированное значение угла), с использованием двух очередных псевдослучайных чисел, полученных с помощью датчика. Таким образом, при переходе к очередному смежному фрагменту выбирали только один (следующий) торсионный угол. Отметим, что в описанном алгоритме номера торсионных углов могут следовать не подряд вдоль по цепи, в общем случае используется последовательность только варьируемых углов.
Если бы при генерировании конформации молекулы в рамках описанной процедуры удалось воспроизвести требуемую плотность
N-з
(1/^зНог^ " ГТ ехр[ ит^ (у7, у7+1, у7+2)/кВТ]
Т=1 __
строго, то искомая оценка Нш, полученная методом МК, имела бы вид, аналогичный выражению (4), т. е.
Ё H (y1, у2 , ..., y5v-1) ■ exP[-Ulong (yV, y2, ..., y5v-1 )/kBT ]
н ш = —-------------s--------------------------------------------------. (10)
Ё eXp[-Ulong (y1, yV , ..., yN-1)/kBT]
V=1
Но состояния в "кубах" вероятностей осуществления углов внутреннего вращения фрагментов задаются лишь с точностью до размеров параллелепипедов. Следовательно, необходимо применять выражение (9), или, в общем виде, выражения (6), (7), в которых следует вычислить реально использованную плотность вероятностей генерирования конформаций р. Это возможно (см. следующий раздел), поскольку размеры параллелепипедов, выбранных при генерировании очередной конформации цепной молекулы, известны. При этом соблюдается и общая рекомендация [12], согласно которой желательно выбирать плотность р по возможности пропорциональной подынтегральной функции, поэтому в итоге можно ожидать небольших дисперсий для оценок Нш.
Плотность вероятностей
Размеры каждого параллелепипеда, который выбирается в процессе генерирования конформации цепи, можно идентифицировать
индексом фрагмента ш7 и тремя случайными номерами Л1)7, Л^ , Лз 7 сторон 1, 2 и 3 конкретного параллелепипеда для каждого номера фрагмента 7 (т. е. торсионного угла) в каждой конформации V данной цепи, где
Л\а = 1, 2, ..., 100; А2)7 = 1, 2, ..., 100; Л^
= 1, 2, ..., 10°. Пусть (^1)то7,Л\п, (£2)тт,Л^Г1, (^з)т7,Л% - длины сторон, в угловых едини-
цах, параллелепипеда, выбранного случайно в конформации номер V во фрагменте номер 7 типа т. Согласно предложенному алгоритму, фиксированные значения углов у1, у2, уз в концевом (т. е. первом, 7 = 1) фрагменте выбраны равновероятно в объеме (¥1)т1д^ ■ (£2)^1 Л21 ■ (рз)т1,Л^,1, поэтому плотность вероятностей генерирования точки с фиксированными значениями углов (у1, у2, уз) равна
100 з ■ [(£1)т1,Л^,1 ■ (£2)т1,Л^,1 ' (£з)т1,Л^,1 ] 1.
Множитель 100-з обусловлен тем, что в данной процедуре выбирается также один из 100
номеров каждой из 3-х сторон параллелепипеда. Генерирование значения угла у4 производится по фрагменту 7 = 2, в котором есть три торсионных угла: у2,Уз,У4, но первые два уже избраны на предыдущем шаге. Плотность вероятностей выбора у4 есть 100-1 ■ [(£з)т2,д£2]-1. Аналогичный сомножи-
тель появляется для всех последующих углов во всех фрагментах вдоль по цепи в конформации с номером V, т. е. для очередного угла номер 7 это 100-1 -[(Ь3)тт,д£ ]-1. Итоговая плотность вероятностей р, с которой реально осуществляется генерирование данной конформации, равна
р = юо-(^-1) • [(і1)ті)Лї! • (і2)ті)Л£! • (¥з)
В общем случае, если варьируются не все углы, в этой формуле в произведении вместо N — 3 должно быть число ^ N — 3. Вы-
ражение р в итоге является простым по сути: в него входят лишь произведения размеров тех параллелепипедов, внутри которых были избраны значения углов. Очевидно, что р из (11) является плотностью, допустимой по отношению к любым функциям, поскольку, в силу конечности размеров всех параллелепипедов, р > 0.
Оценка средних
Для того чтобы осуществить варьирование не только торсионных углов, но и валентных углов и валентных связей (варьирование обобщенное и независимое для каждого валентного угла и связи), после генерирования совокупности торсионных углов проводили, как упоминалось выше, варьирование декартовых коор-
N-3
тх,А5,!]-1 ■ П [(^з)™,А, т]-1. (11)
7=2
динат ж^ ж2,..., жп всех атомов молекулы. При решении ряда специальных задач эту процедуру не осуществляли. Оценка величины диапазона варьирования декартовых координат атомов проведена с использованием параметров поля СнЛИММ27 в варианте [21].
После завершения этапа варьирования декартовых координат всех атомов данной конформации V и получения измененных наборов значений {7}^, {0}^, |у}^ проводили расчет полной энергии молекулы с использованием всех компонентов энергии поля СНЛИММ27 [21] и соответствующих Больцмановских множителей.
Формула (11) дает плотность вероятностей р генерирования конформации, полученной перед проведением этапа варьирования декартовых координат атомов. Именно это выражение следует подставить в (9) или (6), (7) вместо р для получения оценки Нш. Например, используя выражение (7), получим
Е Н* • ехр[-и*/квТ] • [(¥і)ті ^=1
Н ,
Ё ехр[ иу/квТ] ■ [(^1)т1,А? ^=1 ,
Если варьирование декартовых координат
атомов при генерировании конформаций не проводится, то в оценке (12) использовано точное выражение для плотности вероятностей р генерирования конформаций. Эта плотность, хотя и не равна строго требуемой плотности N-3
(1/^«^ог4) 'П вжр[ (у7,у7+1,у7+2)/кВ
7=1
но близка к ней. В итоге выражение (12) да-
— р
ет состоятельную оценку, Нш ——< Н > при ш — те.
Можно продемонстрировать степень эффективности (выгодности) разработанного алгоритма существенной выборки строго по торсионным углам по сравнению с алгоритмом простой выборки, осуществив 2 варианта ((г) и (гг)) компьютерного эксперимента, в ко-
Л^д (р2)ті,Л^,1 (р3)ті,Л3,1 ] • П (^з) тт ,Лз ’ ’ ’ 7=2
N-3
(12)
і (^2)ті,Л^,і (^3)ті,Лз,і ] • П М ,Лз -’ ’ 7=2 ’
торых варьирование декартовых координат атомов не проводится, - валентные связи и валентные углы фиксированы в равновесных положениях.
Ниже приведены оценки ряда средних характеристик для неразветвленной углеводородной цепи без двойных связей, содержащей 10 углеродных атомов (условное обозначение цепи 10:0). Расчеты проведены при температуре Т = 293 К в так называемых 0
- условиях [13] молекулы. Доверительные интервалы во всех приведенных ниже результатах отвечают 95 % надежности согласно ^распределению Стьюдента.
Вариант (і). Цепь 10:0, простая выборка (без варьирования декартовых координат), ш = 9, 6 • 1010 конформаций. Оценки средних значений получены из выражения (9), в ко-
тором плотность р избрана постоянной, случайные значения торсионных углов распределены равномерно: расстояние между концевыми атомами углерода 9,345 ± 0,003 А, квадрат расстояния между концевыми атомами углерода 88,48 ± 0,06 Л2, радиус инерции 3,3434 ±
0,0006 Л.
Вариант (гг). Цепь 10:0, существенная выборка (без варьирования декартовых координат), ш = 7, 2 ■ 108 конформаций. Оценки средних значений получены из выражения (9) с плотностью р распределения торсионных углов в генерируемой выборке согласно (11): расстояние между концевыми атомами углерода 9,344 ± 0,001 Л, квадрат расстояния между концевыми атомами углерода 88,45 ± 0,02 Л2, радиус инерции 3,3433 ± 0,0002 А.
Легко видеть, что оценки всех соответственных характеристик в (г) и (гг) в пределах доверительных интервалов совпадают, но в варианте (гг) использован примерно в 130 раз меньший размер выборки ш, чем в варианте (г), и при этом в 3 раза меньшим оказался доверительный интервал. Если учесть, что вероятная ошибка ~ ш-1/2 [3, 12], то для достижения одинаковых ошибок в (г) и (гг) можно было бы уменьшить ш в варианте (гг) еще примерно в 9 раз. В итоге ш для существенной выборки можно избрать примерно в 1000 раз меньше, чем для выборки простой.
Далее, расчеты показали, что суммарные времена Ь генерирования случайных точек и вычисления всех заданных случайных величин в обоих алгоритмах соизмеримы: на такой расчет для одной конформации цепи по алгоритму (гг) потребовалось всего в 1,25 раза больше времени, чем по алгоритму (г). В итоге произведение времени Ь на ш в варианте существенной выборки (гг) оказалось в несколько сотен раз меньше, чем в варианте простой выборки (г) и, таким образом, разработанный алгоритм существенной выборки является намного более эффективным, выгодным.
Варианты, в которых осуществляется варьирование декартовых координат всех атомов, требуют дополнительного анализа: такое варьирование приводит к изменению не только валентных связей и валентных углов (по отношению к своим равновесным значениям), но и к дополнительному изменению торсионных углов (по отношению к значениям, которые уже были получены в процедуре существенной выборки). Это изменение в формуле
(11) не учтено, и для таких условий формула (11) дает лишь приблизительное выражение для плотности вероятностей р генерирования данной конформации.
В этой связи отметим, что проводимая в настоящей работе последовательность действий может быть рассмотрена и как приложение иных известных способов уменьшения дисперсии [12]. В частности, вычисление интегралов при фиксированных ("замороженных" ) значениях всех величин {1} и {0}, осуществленное на 1-м этапе, фактически означает "интегрирование по части переменных (понижение порядка интеграла)" [12]. Доказано, что если аналитически взять интеграл по некоторым из переменных, а по остальным переменным использовать тот же метод МК, то дисперсия уменьшится. В данном случае этот этап проведен, и можно ожидать, что "размораживание" величин {1} и {0} на 2-м этапе приведет к некоторому увеличению дисперсии.
Далее, варьирование декартовых координат атомов (которое эквивалентно "размораживанию" {1} и {0}), проводится лишь в малом диапазоне (± 0,02 А) вокруг равновесных положений атомов; это отвечает малым колебаниям последних. Ограничение диапазона варьирования декартовых координат до конечных значений означает применение еще одного известного способа уменьшения дисперсии: "интегрирования по части области" [12]. Идея способа такова: если возможно (аналитически) вычислить интеграл по некоторой части области (допустим, что результат интегрирования равен Ш), то всегда выгодно вычислять простым методом МК только интеграл по оставшейся области и добавить Ш к итоговой оценке в качестве слагаемого. Этот способ для данного случая варьирования валентных углов и связей дает Ш — 0, поскольку интегрирование проводится по той части области О, в которой энергия валентных углов и связей — те, т. е. Больцмановский множитель в интеграле — 0. Это есть вся область с большими отклонениями углов и связей. В итоге вычислению методом МК подлежит только интеграл по оставшейся области - вблизи положений равновесия атомов.
Для вычисления интеграла по этой оставшейся области можно рассчитать плотность р, - она будет иметь достаточно сложный вид. Более целесообразно, однако, использовать выражение (11) для генерирования торсионных углов, и рекомендуемый простой метод МК (равномерное распределение) для плотности декартовых координат атомов в узких "потенциальных ямах" (± 0,02 А) вокруг равновесных положений. Типичные значения равновесных длин разных валентных связей в поле СНЛИММ27 [21] -
это ~1,11 - 1,53 А (что отвечает данным физического эксперимента), и амплитуда варьирования ~0,02 А соответствует лишь нескольким процентам от этих величин. Согласно исходным условиям, описанным выше, энергия таких колебаний должна быть сравнительно малой. Естественно ожидать, что влияние таких колебаний и на другие средние характеристики молекул также будет малым.
Методом МК получены оценки степени влияния малых колебаний на ряд характеристик. Для этого осуществлены компьютерные эксперименты (iii) и (iv) с варьированием декартовых координат всех атомов для молекулы 10:0 при T = 293 K в 0 - условиях [13].
Вариант (iii). Цепь 10:0, простая выборка (с варьированием декартовых координат), ш = 9,6 ■ 1010 конформаций. Оценки средних значений получены из выражения (6) (или (7)), в котором положено р = Const, т. е. точки выбирались равномерно в пространстве торсионных углов и в малых областях вокруг равновесных значений атомов: расстояние между концевыми атомами углерода 9,347 ± 0,007 А, квадрат расстояния между концевыми атомами углерода 88,51 ± 0,10 А2, радиус инерции
3,3441 ± 0,0012 A. Эти оценки Нш —>< H > при ш ^ те, поскольку все состояния выбирали равномерно, а дополнительное варьирование торсионных углов при этом не нарушает равномерности распределения по торсионным углам.
Вариант (iv). Цепь 10:0, существенная выборка (с варьированием декартовых координат), ш = 28, 8-10® конформаций. Оценки средних значений получены из выражения (12), т. е. торсионные углы выбирали неравномерно, с плотностью р, но декартовые координаты в пределах областей ± 0,02 А вокруг равновесных значений атомов - равномерно: расстояние между концевыми атомами углерода 9,348 ± 0,001 A, квадрат расстояния между концевыми атомами углерода 88,54 ± 0,02 А2, радиус инерции 3,3443 ± 0,0002 А. Размер выборки ш в варианте (iv) был увеличен в 4 раза по сравнению с таковым в варианте (ii) для того, чтобы добиться совпадения доверительных интервалов в (ii) и (iv).
Как и ожидалось, вследствие учета малых колебаний оценки геометрических характеристик в варианте (iii) увеличились по сравнению с (i) (и, соответственно, в варианте (iv) -увеличились по сравнению с (ii)). Вместе с тем относительный рост этих характеристик составил всего 0,02-0,04 %; т. е. влияние обсуждаемого варьирования декартовых координат
атомов данной цепной молекулы на ее геометрические характеристики, действительно, оказалось малым.
Отметим, что оценки всех соответственных характеристик вариантов (ггг) и (г-и) в пределах доверительных интервалов совпадают (аналогично тому, как совпадают оценки характеристик в вариантах (г) и (гг)). Размер выборки ш в варианте (ггг) примерно в 30 раз больше, чем в варианте (г-и), а доверительные интервалы - больше примерно в 5-7 раз, т. е. для получения одинаковых интервалов в (ггг) и (г-и) следует увеличить ш в варианте (ггг) еще примерно в 25-50 раз, и получим в ^750-1500 раз больше, чем в варианте (г-и). Поскольку времена расчета случайных величин соизмеримы, вариант существенной выборки с варьированием координат (г-и) оказался также в сотни раз эффективнее, чем соответствующий вариант простой выборки (ггг).
Разработанный алгоритм компьютерной имитации цепных молекул позволяет воспроизвести точное химическое строение цепи, учитывает взаимозависимость каждых трех последовательных углов внутреннего вращения и непрерывный спектр конформаций цепей, а при генерировании этих конформаций использует метод существенной выборки. В итоге он позволяет передать основные особенности молекул цепного строения и, обладая общностью, может быть адаптирован для конструирования и исследования свойств цепных молекул разных классов.
Работа выполнена при поддержке РФФИ (проект 10-03-00201а) и программы Президента РФ - Ведущие научные школы (НШ-1642.2012.4).
ЛИТЕРАТУРА
1. Биндер К. Методы Монте-Карло в статистической физике / Ред. К. Биндер. М.: Мир, 1982. С. 9-57.
2. Биндер К., Хеерман Д. В. Моделирование методом Монте-Карло в статистической физике: Введение. М.: Наука; Физматлит, 1995. 144 с.
3. Бусленко Н. П., Голенко Д. И., Соболь И. М., и др. Метод статистических испытаний (метод Монте-Карло). М.: ГИФМЛ, 1962. 332 с.
4. Гнеденко Б. В. Курс теории вероятностей. М.: Гос. изд. физ.-мат. лит., 1961. 408 с.
5. Гулд Х., Тобочник Я. Компьютерное моделирование в физике. М.: Мир, 1990. Ч. 2. 400 с.
6. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука; Гл. ред. физ.-мат. лит., 1971. 328 с.
7. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М.: Наука; Гл. ред. физ.-мат. лит., 1982. 296 с.
8. Зализняк В. Е. Основы научных вычислений. Введение в численные методы для физиков. М.: Едиториал УРРС, 2002. 296 с.
9. Ландау Л. Д., Лифшиц Е. М. Статистическая физика. Ч. 1. М.: Наука; Гл. ред. физ.-мат. лит., 1976. 584 с.
10. Методы компьютерного моделирования для исследования полимеров и биополимеров / Отв. ред. В. А. Иванов, А. Л. Рабинович, А. Р. Хохлов. М.: Книжный дом ЛИБРОКОМ, 2009. 696 с.
11. Рабинович А. Л., Иванов В. А. Методы компьютерного моделирования для исследования полимеров и биополимеров / Отв. ред. В. А. Иванов, А. Л. Рабинович, А. Р. Хохлов. М.: Книжный дом ЛИБРОКОМ, 2009. С. 63-119.
12. Соболь И. М. Численные методы Монте-Карло. М.: Наука; Гл. ред. физ.-мат. лит., 1973. 312 с.
13. Флори П. Статистическая механика цепных молекул. М.: Мир, 1971. 440 с.
14. Хеерман Д. В. Методы компьютерного эксперимента в теоретической физике. М.: Наука; Гл. ред. физ.-мат. лит., 1990. 176 с.
15. Allen M. P., Tildesley D. J. Computer Simulation of Liquids. Oxford: Clarendon Press, 1987. 385 p.
16. Feller S. E., MacKerell, Jr. A. D. An improved empirical potential energy function for molecular simulations of phospholipids // J. Phys. Chem. B. 2000. Vol. 104. P. 7510-7515.
17. Feller S. E., Yin D., Pastor R. W, MacKerell, Jr. A. D. Molecular dynamics simulation of unsaturated lipids at low hydration: parametrization and comparison with diffraction studies // Biophys. J. 1997. Vol. 73. P. 2269-2279.
18. Flory P. J. Foundations of rotational isomeric state theory and general methods for generating
СВЕДЕНИЯ ОБ АВТОРАХ:
Рабинович Александр Львович
гл. науч. сотр., д. ф.-м. н.
Институт биологии Карельского научного центра РЛИ ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: rabinov@krc.karelia.ru тел.: (8142) 571879
Ж^уркин Дмитрий Викторович
аспирант
Петрозаводский государственный университет, физико-технический факультет
ул. Университетская, 10 А, Петрозаводск, Республика
Карелия, Россия, 185910
эл. почта: zhurkin_dm@mail.ru
тел.: (8142) 719671
configurational averages // Macromolecules. 1974. Vol. 7, N 3. P. 381-392.
19. Frenkel D., Smit B. Understanding Molecular Simulation. From Algorithms to Applications. 2nd edition, San Diego: Academic Press, 2002. 638 p.
20. Go N., Scheraga H. A. On the use of classical statistical mechanics in the treatment of polymer chain conformations // Macromolecules. 1976. Vol.
9, N 4. P. 535-542.
21. Hogberg C.-J., Nikitin A. M., Lyubartsev A. P. Modification of the CHARMM force field for DMPC lipid bilayer // J. Comput. Chem. 2008. Vol. 29. P. 2359-2369.
22. Kahn H. Random sampling (Monte Carlo) techniques in neutron attenuation problems. 1 // Nucleonics. 1950. Vol. 6, N 5. P. 27-33.
23. Kahn H. Random sampling (Monte Carlo) techniques in neutron attenuation problems. 2 // Nucleonics. 1950. Vol. 6, N 6. P. 60-65.
24. Klauda J. B., Brooks B. R., MacKerell, Jr. A. D. et al. An ab initio study on the torsional surface of alkanes and its effect on molecular simulations of alkanes and DPPC bilayer // J. Phys. Chem. B. 2005. Vol. 109. P. 5300-5311.
25. Landau D. P., Binder K. A Guide to Monte Carlo Simulations in Statistical Physics. 2nd edition. Cambridge University Press, 2005. 432 c.
26. Simulation Methods for Polymers / Eds. Kotelyanskii M., Theodorou D.N. N.Y.: Marcel Decker, Inc., 2004. 602 p.
27. Matsumoto M., Nishimura T. Mersenne
Twister: A 623-dimensionally equidistributed
uniform pseudorandom number generator // ACM Transactions on Modeling and Computer Simulations: Special Issue on Uniform Random Number Generation. 1998. Vol. 8. P. 1-25.
28. Schlenkrich M., Brickmann J., MacKerell, Jr. A. D., Karplus M. In: Biological membranes: a molecular perspective from computation and experiment / K. M. Merz, B. Roux, Eds. Birkhauser, Boston. 1996. P. 31-81.
Rabinovich, Alexandr
Institute of Biology, Karelian Research Centre, Russian Academy of Sciences
11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: rabinov@krc.karelia.ru tel.: (8142) 571879
Zhurkin, Dmitrii
Petrozavodsk State University, Faculty of Physical Engineering
10 A University St., 185910 Petrozavodsk, Karelia, Russia e-mail: zhurkin_dm@mail.ru tel.: (8142) 719671
