Моделирование н-алканов и полиглицинов методом Монте-Карло с использованием алгоритма Ванга-Ландау Текст научной статьи по специальности «Физика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 4 2007 Вып. 3

УДК 537.213

А. А. Юрченко, П. Н. Воронцов-Вельяминов

МОДЕЛИРОВАНИЕ Н-АЛКАНОВ И ПОЛИГЛИЦИНОВ МЕТОДОМ МОНТЕ-КАРЛО С ИСПОЛЬЗОВАНИЕМ АЛГОРИТМА ВАНГА-ЛАНДАУ

Введение. Метод Монте-Карло (МК), предложенный пятьдесят лет назад Метро-полисом и др. [1], оказался мощным инструментом в изучении неидеальных молекулярных систем [2, 3]. Однако стандартная МК-процедура становится неэффективной во множестве важных физических ситуаций. Для изучения систем, имеющих грубый рельеф потенциальной энергии с множеством локальных минимумов, при рассмотрении фазовых переходов и других явлений, проходящих при низких температурах, высоких плотностях или в присутствии сложных молекулярных компонентов, приходится модифицировать стандартные подходы. Такие модификации, предложенные в прошлом десятилетии, известны теперь под общим названием обобщенных ансамблей МК [4, 5]. К ним относятся расширенный ансамбль МК [6, 7] и энтропическое моделирование [8, 9]. Эти методы оказались эффективными в решении перечисленных выше проблем, хотя и у них имеется общий недостаток: они требуют предварительной настройки набора параметров («балансирующих факторов») [6], играющих ключевую роль при моделировании. В 2001 г. Ванг и Ландау (ВЛ) предложили алгоритм, в котором происходит автоматическая настройка параметров энтропического моделирования [10].

В наших предыдущих работах [11, 12] мы применили ВЛ-алгоритм к энтропичес-кому моделированию простых полимерных моделей: свободносочлененных свободных и замкнутых цепей в трехмерном континуальном пространстве и свободных и замкнутых цепей на простой кубической решетке. В атермическом случае для континуальных моделей нами был предложен вариант ВЛ-алгоритма, позволивший определить избыточную энтропию полимера как функцию диаметра мономера во всем диапазоне его изменения в одной МК-процедуре. В термическом случае были вычислены распределения по энергиям (для решеточных полимеров—по числу контактов). Эти распределения позволили рассчитать значения конформационной внутренней энергии, теплоемкости, избыточной энтропии и среднеквадратичные значения радиуса инерции. В настоящей работе ранее развитые алгоритмы применяются для изучения моделей цепи с фиксированным валентным углом—предельных углеводородов и полиглицина.

1. Используемые модели и алгоритмы

1.1. Метод Монте-Карло в классической статистической механике. Задачи равновесной статистической механики классических систем можно свести к вычислению статистического интеграла в каноническом ансамбле [13], который после интегрирования по импульсам приводится к виду Z(N,V,T) = -~j^j"exp{-ߣ/(<7)}i&7, где Л = -y/^f-— тепловая длина волны де Бройля частиц массы т при температуре Т. Таким образом, задача сводится к вычислению конфигурационного интеграла ß(ß) = Jexp{-ß{/(<?)}£&?.

О A.A. Юрченко, П.Н. Воронцов-Вельяминов, 2007

От интегрирования по координатам можно перейти к интегрированию по энергии:

- Е)йц, где 0.(Е)—объем части конфигурационного пространства, в которой потенциальная энергия системы лежит в пределах от £ до Е + с1Е. Вычисление С1(Е) в общем случае возможно только численными методами. Диапазон изменения энергии системы ЕтЫ <Е< Етах разбивается на конечное число (Лу равных отрезков («ящиков»), и в вычислительном эксперименте определяются доли £1{Е), соответствующие г'-му отрезку энергии. При вычислении канонических средних от интегралов переходят к интегральным суммам. Среднее каноническое от вели-

Nь ¡Нь

чины /получают по следующей формуле (/)(Р) = ехР{~Р^} / А ехР(~Р£,}>

¿=1 / 1=1

где ^— значение величины / для /-го отрезка энергии.

1.2. Алгоритм Ванга-Ландау. Алгоритм Ванга-Ландау [10] является реализацией алгоритма энтропического моделирования [9]. Он решает проблему подбора подходящих вероятностей перехода для получения требуемого при энтропическом моделировании равномерного посещения энергетических состояний и, следовательно, позволяет получить плотность состояний 0.(Е),

Диапазон изменения энергии системы Е <Е< Етах разбивается на конечное число (Лу равных отрезков («ящиков»). Заводится массив Г2, состоящий из N элементов, каждый из которых соответствует отрезку разбиения энергии. Изначально все элементы Ц берутся равными единице. В процессе вычислительного эксперимента на каждом МК-шаге происходит изменение конформации системы. Пусть Ех и Ег—это энергии системы до изменения и после. Каждая из них попадает в свой «ящик»—/-й и у'-й, соответственно (номера / и ] могут совпадать). В таком случае, изменения в системе принимаются с вероятностью

р(Е{ -> Е2) = тш

в случае отказа система возвращается в исходное состояние. После принятия или не принятия новой конформации системы все повторяется на новом МК-шаге. Каждый раз при посещении к-го «ящика» (в случае принятия изменений системы к = ], при отказе к = /) производится изменение к-го элемента массива П. Он умножается на инкремент а > 1, т. е. —► • а. На протяжении серии МК-шагов величина инкремента остается неизменной. На каждой последующей серии значение параметра а уменьшается. В работе Ванга и Ландау [10] использовалось рекуррентное соотношение а0 = е ~ 2,71828, ат+1 - где т—номер серии. Отметим, что для модификации величины а подходит любая функция, которая монотонно стремится к единице. В результате этого алгоритма происходит автоматическая настройка весов вероятности перехода (*), которые одновременно являются плотностями состояний. По окончании вычислительного экс-

перимента производится нормирование массива Г2 на единицу: 0.к Полу-

чившийся массив П является функцией распределения по энергиям.

Одновременно с массивом О. заводится массив посещений V, элементы которого изначально равны нулю. На каждом МК-шаге в ячейку Ук, соответствующую по-

сещению к-го «ящика», добавляется единица. В каждой серии добиваются равномерной гистограммы посещений с достаточной степенью точности.

Если система может принимать состояния с энергией меньше Е или больше Етах, заводятся дополнительные «ящики»: один для Е<ЕтЫ, второй для Е>Етах. При наличии самопересечений системы, этим состояниям ставится в соответствие еще один «ящик», а если доля самопересекающихся состояний велика, то несколько «ящиков», соответствующих разному числу самопересечений. Алгоритм при этом не меняется: на каждом МК-шаге система «переходит» из одного «ящика» в другой, этот переход принимается с вероятностью (*), изменяются соответствующие элементы массивов О и V, и все повторяется на новом МК-шаге.

Описанную выше процедуру можно сделать более удобной, если перейти к энтропии состояний 8(Е) = 1пГ2(£'). Вместо массива Ф заводится массив энтропий Я. В таком случае вероятность перехода (*) перепишется как р(Е1 —» Е2) = тт(1,ехр{5, На каждом МК-шаге энтропия принятого состояния будет изменяться Бк —* Бк + Д5т, где = 1п ат, т—номер серии. По окончании вычислительного эксперимента функция

распределения по энергиям вычисляется по формуле: 0.к = ехр(^) / ^Г ехр(5/). В насто-

/

ящей работе величина Д50 полагалась равной 0,1. А в каждой следующей серии значение Д5 уменьшалось по закону Д5т = Д50 ■ 0,75™, где т—номер серии. Всего проводилось 25 серий, в которых число МК-шагов постепенно увеличивалось, и всего проводилось около 20 млн МК-шагов. Относительное отклонение гистограммы посещений от равномерной не превышало 10%.

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

1.3. Модель предельных углеводородов. В первой части работы исследованы цепи с фиксированным валентным углом—модель неразветвленных предельных углеводородов СН^СНу^-СН,. Цепь состоит из И+ 1 мономеров, соединенных N сегментами длиной /=0,1526 нм. На концах цепи расположены мономеры СН^ во всех остальных узлах—мономеры СН2. Все пары соседних сегментов образуют между собой фиксированный угол 109°28\ В цепи разрешено свободное вращение вокруг любого сегмента [14].

Построение цепи осуществлялось следующим образом. Генерировалось Ы- 1 случайных переменных ц>к(к=2^ 3,... №), равномерно распределенных на отрезке [0,2л]. Каждая из них в паре с углом 0^= 70°32' являлись локальными сферическими координатами центра к-го мономера в системе координат, в которой центр (к- 1)-го мономера находится в начале координат, а центр (к-2)-го мономера—в точке (0,0,-1). Далее все локальные координаты переводились в лабораторные декартовые координаты. В лабораторной системе отсчета центр нулевого мономера располагали в начале координат, центр первого мономера—в точке (0,0, Г). Чтобы определить координаты остальных мономеров, записывали координаты каждого мономера относительно предыдущего Хк= (0,0, Г) и затем переводили их в лабораторные координаты, используя обратную матрицу преобразования А~1 [15] и сдвиг вдоль оси г:Х0 =_(0, 0, 0); Х1 = Х;,Х2 = А?Х2 + X,; Х3 = А? (А;хХг +Х1) + Х1 ... Хк = А?(...(Ак^(Ак'Xк + Хк_х) + Хк_г)...) + Хх. Конформацию цепи изменяли поворотом части цепи, находящейся от случайно выбранного сегмента до ее конца, вокруг оси, которая проходит через выбранный сегмент на случайный угол, выбираемый равномерно от 0 до 2л. В результате все валентные углы оставались без изменения.

Конформации полимера сортировались по величине потенциальной энергии

V = ^и,ог + где ^и^—суммарная энергия торсионных взаимодействий между всеми парами мономеров, расположенных через два мономера, —суммарная энергия взаимодействий между всеми парами мономеров, расположенными через три и более мономеров. Энергия торсионного взаимодействия определялась по формуле и,ог =-10,66-16,79со8ф + 1,13со82ф + 26,32со83ф (кДж/молъ) [14], где ф—угол между плоскостями, проходящими через пары сегментов, соединяющих взаимодействующие мономеры. Энергия взаимодействия между мономерами, расположенных через три и более мономеров, определялась потенциалом Леннард—Джонса с бесконечным «барьером» (бесконечная энергия соответствовала самопересечению цепи):

и» =

4е

( 12 _ л

\ и и У

г>ст

оо, г <ст

где г—расстояние между центрами мономеров, а = 0,3905 нм, е/у = ф п, величина е., для мономера СН3 равна 0,733 кДж/моль, для мономера СН2—0,494 кДж/моль [14].

Диапазон изменения С/ определялся в предварительном расчете. Этот отрезок энергий ЕтЫ <и< Етах разбивался на Иь = 100 «ящиков». В результате ВЛ-вычисли-тельного эксперимента определялась функция распределения по «ящикам» (О.) и далее были рассчитаны значения внутренней энергии, теплоемкости, радиуса инерции и квадрата расстояния между концами цепи для различных длин полимера (5 <N<30).

1.4. Модель полипептидов. Полипептидные цепи моделировались на атомарном уровне—молекула определялась как совокупность атомов водорода (Н), углерода (С), кислорода (О) и азота (Ъ1). Расстояния и углы между ними задавались из базы данных СНАИММ [16]. Длина валентных связей и величина валентных углов фиксировались. В качестве радикалов при всех Са брались водороды -Н, т. е. исследовался полиглицин. Таким образом, в качестве модели использовалась цепь с фиксированными валентными углами, в которой разрешались вращений только вокруг связей 1Ч-Са и Са-С.

Построение начальной конформации полимера осуществлялось схожим образом с построением цепи с валентным углом (п. 1.3), однако, с учетом того, что расстояния между различными типами атомов отличаются друг от друга. Вначале происходило построение полипептидного остова (-1Ч-Са-С-)и, а затем—достраивание атомов О и Н к атомам углерода и азота. Изменение конформации цепи осуществляли поворотом части молекулы вокруг оси, проходящей через случайно выбранную связь (вокруг которой разрешено вращение), на случайный угол, выбираемый равномерно от 0 до 2я. Конформации цепей сортировали по величине потенциальной энергии и, которую определяли как сумму потенциалов Леннард-Джонса между всеми парами атомов, находящихся на нефиксированном расстоянии. Константы взаимодействия брались из базы данных СНАКММ [16]. ВЛ-вычислительный эксперимент проводился для различных длин полимера (2<п< 10).

2. Результаты и их обсуждение

2.1. Модель предельных углеводородов. Изучалась модель предельных углеводородов, параметры которой описаны в п. 1.3. В результате были определены функции распределения по энергиям для полимеров различной длины (рис. 1). Под

Еуд (кДж /моль )

Рис. 1. Нормированные функции распределения по удельным энергиям модели предельных углеводородов для полимеров различной длины Ы= 6 (♦), 9 (А), 12 (■), 16 (•), 18 (0), 19 (А), 20 (□), 24 (о) и 30 (х) сегментов/

удельной энергией понимается полная энергия, поделенная на N—2—число торсионных взаимодействий, которое совпадает с числом степеней свободы системы. Общие доли самонепересекающихся цепей С1БА}У в зависимости от длины цепи представлены в табл. 1. Поскольку эти доли очень малы (особенно для длинных полимеров), а доли самопересекающихся цепей велики, производилась сортировка по числу самопересечений. Заводилось четыре «ящика», первый из которых соответствовал цепям с 1-2 самопересечениями, второй—3-5, третий—6-10, четвертый—с более чем 10 самопересечениями. Это позволило избежать «зависания» системы в одном «ящике», соответствующем много большему объему конформационного пространства, чем остальные.

Полученные функции распределения Ц позволили рассчитать значения внутренней энергии и теплоемкости, средние квадраты радиусов инерции и расстояний между концами полимера. На рис. 2 представлены зависимости средней энергии от температуры (сплошные линии). Видно, что для всех длин полимеров средняя каноническая энергия монотонно возрастает с увеличением температуры. Открытыми квадратиками нанесены значения Т(Е), полученные следующим образом. Известно, что Т - (§§-)^, [13], а 5 = И\п£1. Тогда, переходя к конечным разностям, из зависимостей

Таблица 1

Зависимость доли самонепересекающихся цепей модели предельных углеводородов от числа сегментов (Л/)

N 6 9 12 16 18 19 20 24 30

^АЯ' 2,3-Ю-1 5,8-Ю-2 1,7-10"2 2,7-10 3 1,0-ю-3 6,2-10"4 3,9-Ю^1 6,6-ю-5 4,1-10"®

Рис. 2. Зависимости внутренней энергии модели предельных углеводородов от температуры для полимеров различной длины N = 6, 9, 12, 16, 18, 19,20, 24 и 30 сегментов.

Открытыми квадратиками нанесены значения функций Т(Е) = д полученных из представленных на рис. 1 данных.

£1(Е) можно определить температуру как функцию от энергии: Т(Е) = ^—¡у. Функции энергии от температуры, полученные независимыми друг от друга интегральным и дифференциальным методами, прекрасно согласуются между собой (рис. 2).

На зависимостях теплоемкости от температуры наблюдается стремление С(Т) к нулю при Г—* 0 и Г—► оо. Для цепей длиной N<20 наблюдается совпадение положений максимумов—все они находятся в пределах 100 ±20 К. Для более длинных цепей происходит сдвиг максимума в большие температуры.

Для полимеров длиной Ы> 20 наблюдается увеличение среднеквадратичного радиуса инерции при росте температуры, что может свидетельствовать о переходе «клубок—глобула». В то же время для коротких полимеров (7У< 20) наблюдается обратная картина: увеличение температуры приводит к уменьшению среднеквадратичного радиуса инерции. Это связано с фактом наибольшего вклада в энергию торсионного взаимодействия. При низких температурах полимер стремится приобрести конформацию с наименьшей энергией.

Таблица 2 Сравнение среднеквадратичных расстояний между концами полимеров (<Иг>ш) с определяемым из геометрических соображений расстоянием (//) между концами цепочки, все связи которой находятся в транс-состоянии, для полимеров с различным числом сегментов (А)

N И (А) <А2>Ш (А)

6 1,476 7,475

9 11,25 11,23

12 14,95 14,93

16 19,94 19,89

18 22,43 22,36

19 23,69 20,89

20 24,92 12,55

24 29,90 9,44

30 37,38 12,29

Рис. 3. Конформации модели предельных углеводородов длиной N= 12 (а, д, и), 18 (б, е, к), 19 (в, ж, л) и 24 (г, з, м) сегментов с различными энергиями: -65 (а), -130 (д), -195 (и), -100 (б), -185 (е), -315 (к), -110 (в), -195 (ж), -325 (л), -135 (г), -270 (з) и -435 кДж/моль (м).

И при превалировании вклада торсионного взаимодействия этой конформацией становится вытянутая цепочка, все связи которой лежат в одной плоскости и находятся в транс положении. Действительно, если сравнить среднеквадратичные расстояния между концами полимеров при стремлении темпераггуры к нулю с расстояниями между концами в транс-цепочках (табл. 2), то можно обнаружить совпадение этих величин для всех длин полимеров N<18, небольшое различие (~ 12%) при Лг= 19 и существенное расхождение при N>20. Таким образом, можно придти к выводу, что при низких температурах короткие полимеры вытягиваются в палочкообразную структуру, а более длинные полимеры образуют компактные структуры.

На рис.3 представлены конформации полимеров длиной N=12, 18, 19 и 24 сегмента при различных энергиях. Можно видеть, что цепи длиной Ы = 12, 18 и 19 сегментов с самыми низкими энергиями, а именно эти конформации стремятся принять полимеры при уменьшении температуры к нулю, имеют вытянутую палочкообразную структуру. В то же время полимер длиной N= 24 сегмента в низкоэнергетической конформации, которую он достигает в МК-процедуре, образует компактную структуру, которая имеет палочкообразные части.

Для высоких температур наблюдается монотонный рост среднеквадратичных радиусов инерции и расстояния между концами полимера при увеличении его длины.

Таблица 3

Зависимость предельного значения среднеквадратичного расстояния между концами полимера (<й2>1/2) при Тоо от числа сегментов (АО

N 6 9 12 16 18 19 20 24 30

<иг>ш (А) 6,27 8,88 11,23 14,05 15,33 15,98 16,58 18,85 22,00

В табл. 3 представлены предельные значения среднеквадратичного расстояния между концами цепи при Г —> со. При больших температурах полимер не является фантомным в силу запрета на самопересечение. В связи с этим зависимость <к2>т{И) не является линейной. Приведенные в табл. 3 данные ложатся на степенную функцию от .¡V с показателем степени 0,77 и квадратом смешанной корреляции Я2 = 0,99822. В то же время, если рассматривать зависимость <И2>т(Ы- 2), то данные укладываются на степенную функцию с показателем степени 0,65 и квадратом смешанной корреляции Я2 = 0,99987.

2.2. Модель полипептидов. Исследовалась модель полиглицина, параметры которой описаны в п. 1.4. Рассчитаны функции распределения по энергиям для числа мономеров-глицинов п = 2, ,.., 10. Общие доли самонепересекающихся полимеров 0.ЗАЦ, в зависимости от длины цепи представлены в табл. 4. Полученные функции распределения П позволили рассчитать значения внутренней энергии, теплоемкости и средние квадраты радиусов инерции полимеров. На рис. 4 представлены зависимости средней энергии от температуры. Видно, что для всех п средняя каноническая энергия монотонно возрастает с увеличением температуры. При низких температурах (Г—>0) наблюдается монотонное уменьшение средней канонической энергии с увеличением

Т(К)

Рис. 4. Зависимости внутренней энергии модели полиглицина от температуры для полимеров длиной от 2 до 10 пептидов. Серым цветом выделена область от 250 до 350 К.

Рис. 5. Конформации полиглицина с различным числом п- 5 (а-в) и 10 (г-д) мономеров-глицинов с различными энергиями: 126 (а), 73 (б), -31 (в), 126 (г), 42 (д) и -84 (е) кДж/моль.

числа мономеров, в то время как для высоких температур (Т —*■ со) наблюдается обратная картина: чем больше число мономеров в полиглицине, тем больше средняя каноническая энергия.

С увеличением длины полимера происходит монотонный сдвиг максимума теплоемкости в сторону увеличения температуры. Для всех длин полимеров наблюдается рост среднеквадратичного радиуса инерции при увеличении температуры, что может свидетельствовать о переходе «клубок-глобула». Причем, чем больше число мономеров в полиглицине, тем сильнее выражен этот переход, а его положение сдвигается в большие температуры. Действительно, если построить изображения конформаций полиглицинов при различных энергиях (рис. 5), можно отметить, что с понижением энергии происходит компактизация полимерного клубка.

Заключение. Использованный в работе алгоритм Ванга-Ландау показал свою эффективность. По сравнению со стандартным алгоритмом Метрополиса [1] ВЛ-алго-ритм [10] обладает рядом достоинств. При числе шагов около 107 ВЛ-алгоритм позволяет вычислять доли, значения которых различаются на 20-25 порядков (рис. 1). Функции распределения, полученные в широком диапазоне своих значений, позволяют получать канонические средние в широком диапазоне температур. При этом функция распределения рассчитывается за один вычислительный эксперимент.

В работе исследовались цепи сравнительно малой длины. Для модели предельных углеводородов длиной 30 сегментов и модели полиглицина, состоящей из 10 пептидов,

время вычислительного эксперимента на компьютере с частотой процессора 2 ГГц составило 21 и 75 часов соответственно. Увеличение длины цепи вдвое приведет к увеличению времени счета вчетверо. На обсчет системы потребуются недели счетного времени. Для увеличения производительности можно разбить исследуемый диапазон энергий на несколько отрезков, на каждом из которых с помощью алгоритма Ванга-Ландау рассчитать функцию распределения по энергиям, а затем «сшить» ее на всех отрезках. Уменьшение исследуемого диапазона уменьшит время счета, а использование многих процессоров позволит параллельно рассчитывать систему на выбранных отрезках энергий.

В работе исследовались полипептиды, пептиды которых были выбраны самой простой структуры: радикалы у атомов Сц являлись атомами водорода. Заменяя атом водорода другими радикалами, можно строить модели полипептидов с любой последовательностью пептидов и изучить их конформационные свойства.

Работа выполнена при содействии гранта РФФИ 05-02-17428.

Summary

Yurchenko A. A., Vorontsov-Velyaminov Р. N. Computer simulation of n-alkanes and polyglycines using Monte-Carlo method within Wang-Landau algorythm.

In this article we present computer simulation of polymer chains with fixed valence angle. Monte-Carlo methods within Wang-Landau algorithm was used for simulation of n-alkanes and polyglycines. The obtained distribution functions provided of calculation of energy, heat capacity, radius of inertia and end-to-end distance in a wide temperature range. Forcefield OPLS and CHARMM were used in simulations of n-alkanes and polyglycines respectively,

Литература

1. Metropolis N.. Rosenbluth A. W., Rosenbluth M.N., Teller A. H.,Teller E. И J. Chem. Phys. 1953. Vol. 21. P. 1087-1092. 2. Binder K. Monte Carlo Methods in Statistical Physics. Berlin, 1979. 3. Allen M. P., Tildesley D.J. Computer Simulation of Liquids. Oxford, 1987. 4. IbaY.II Int. J. Mod. Phys. C. 2001. Vol. 12. P. 623. 5. MitsutakeA., Sugita Y., OkamotoY.II Biopolymers. 2001. Vol. 60. P. 96. 6. Lyubartsev A.P., MartsinovskiiA.A., Shevkunov S.V., Vorontsov-Velyaminov P.N. I I J. Chem. Phys. 1992. Vol. 96. P. 1776. 7. MarinariE., Parisi G. //Europhys. Lett. 1992. Vol. 19. P. 451. 8. BergB.A, Neuhaus Т. II Phys. Rev. Lett. 1992. Vol.68. P. 9. 9. Lee J. II Phys. Rev. Lett. 1993. Vol. 71. P. 211. 10. WangF., Landau D. P. II Phys. Rev. Lett. 2001. Vol. 86. P. 2050. 11. Vorontsov-Velyaminov P.N., VolkovN.A., Yurchenko А. А. И J. Phys. A: Math. Gen. 2004. Vol. 37. P. 1573-1588. 12. VolkovN.A., Yurchenko A. A., Lyubartsev A. P., Vorontsov-Velyaminov P.N. II Macromol. Theory Simul. 2005. Vol. 14. P. 491-504. 13. Ландау JI.Д., ЛифшицЕ.М. Статистическая физика. М., 1964. 14. Jorgensen W.L., Madura J.D., Swenson С J. IП. Am. Chem. Soc. 1984. Vol. 106. P. 6638-6646. 15. Goldstein H. Classical Mechanics. Cambridge, 1950. 16. Brooks C.L., Bruccoleri R.E., OlafsonB.D., States D. J., Swamiathas S., Karplus. M. //J. Comput. Chem. 1983. Vol. 4. N 1. P. 187-217.

Статья принята к печати 27 февраля 2007 г.