Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 20
Текст из файла (страница 20)
Моделирование методом классических траекторий требует большого числа разнообразных эффективных процедур вычислений. Совокупность этих процедур представляет собой математическое обеспечение метода классических траекторий. В следующем разделе будут рассмотрены наи. более важные процедуры, используемые при динамических расчетах. 4.
ЧИСЛЕННЫЕ МЕТОДЫ, ИСПОЛЬЗУЕМЫЕ ПРИ ДИНАМИЧЕСКИХ РАСЧЕТАХ Ниже будут проанализированы различные схемы интегрирования классических уравнений движения, рассмотрен метод вычисления пучка близ. ких траекторий и метод глобального исследования поверхностей потенциальной энергии (поиск точек перевала), приведены статистические процедуры, используемые при расчетах методом классических траекторий.
Численное интегрирование систем уравнений Гамильтона При динамических расчетах в тех случаях, когда требуется рассчитывать не слишком длинные траектории (менее 100 колебаний взаимодействующих молекул(, как правило, используются обычные разностные методы решения обыкновенных дифференциальных уравнений. Однако для получения достоверных результатов приходится осуществ- лять интегрирование нескольких тысяч траекторий, и, таким образом, выбранный метод должен обладать максимальным быстродействием. Для выбора наиболее эффективного метода интегрирования системы дифференциальных уравнений Гамильтона необходимо установить: метод должен быть одношаговым (методы Рунге — Кутта) или многошаговым (методы типа предиктор-корректор), с автоматическим выбором шага или с постоянным шагом интегрирования, кроме этого, необходимо найти оптимальный порядок метода.
В ранних работах, выполненных методом классических траекторий, очень популярным был (иетод Рунге-Кутта-Гилла 4 го порядка [265). Гюзже стали применяться многошаговые методы высокого (до 15го) порядка точности, в основном использовалась процедура Адамса-Мултона 4-го порядка [92), метод экстраполяции [219).
В большинстве работ численное решение систем дифференциальных уравнений осуществлялось с постоянным шагом интегрирования. В [66) было проведено сравнение эффективности различных процедур численного интегрирования для решения систем обыкновенных дифференциальных уравнений. Исходя из уже имеющегося опыта решения различных задач метода классических траекторий, можно сделать следующие заключения относительно выбора того или иного способа интегрирования. 1.
Решения системы уравнений Гамильтона являются гладкими функциями, как правило, они имеют осциллирующий характер. С учетом этого целесообразно применять методы высокого порядка точности. С другой стороны, из-за ограниченной точности машинного представления чисел (для ЭВМ БЭСМ 6, например, мантисса числа занимает 40 бит) неовлв сообразно применять методы выше 6 го порядка точности [66, Г74) . 2. Следует использовать методы с автоматическим выбором шага интегрирования. Это обусловлено тем, что в этом случае интегрирование осуществляется с оптимальным шагом, который, естественно, сильно различается для случаев наибольшего сближения атомов (в середине траектории) и наибольшего их разлета (на концах траектории).
Некото. рое усложнение программирования процедуры интегрирования не имеет особого значения, так как исследователь, как правило„пользуется уже созданными эффективными и проверенными пакетами программ. Таким образом, для решения задач методом классических траекторий (и других систем обыкновенных диффцзенциальных уравнений, не обладающих большой жесткостью) можно рекомендовать следующие процедуры и пакеты программ, реализующие разностные методы решения систем обыкновенных дифференциальных уравнений: 1) процедуру КОТТАМЕВЗО(ч' [324), реализующую одношаговый метод Кутта — Мерсона с автоматическим выбором шага, 4.го порядка точности с пятью вычислениями правых частей на каждом шаге интегрирования и контролем локальной величины ошибки по Мерсону; 2) пакет программ ВКР45 [182), реализующий одношаговый метод Рунге-Кутта 4-го порядка, предложенный Фелбергом [259), требующий шести вычислений правых частей на каждом шаге интегрирования.
Как и в процедуре КОТТАМЕВВО(ч', сравнение решений по различным методам 4.го и Ьго порядков используется для контроля точности вьг числений; 3) алгоритм, предложенный в работах [125„343), реализующий многсь шаговый метод Адамса 4-го порядка с автоматическим выбором шага интегрирования; 4) пакет программ ВТЕР и ОЕ [399), реализующий многошаговый метод Адамса с переменным порядком и шагом интегрирования: 77 б) многофункциональнью пакеты программ ЯТ(РЕ [67) и ЕР(ВООЕ [2261, реализующие многошаговые методы Адамса и Гира с переменным порядком и шагом интегрирования. Процедуре КОТТАМЕВЯО(ч наиболее компактна и проста в обраще нии. Согласно результатам работ [66, 1091, где проводилось сравнение различных численных матодов решения систем обыкновенных дифференциальных уравнений, она очень эффективна и ее можно усиленно рекомендовать при расчетах методом классических траекторий.
При моделировании некоторых элементарных процессов (мономолеку. лярный распад многоатомных молекул, процессы пе(жраспределения энергии) возникает необходимость расчета длинных (более 100 низко. частотных колебаний молекульб траекторий. В этом случае оказывается, что применение ревностных методов может привести к существенному накоплению ошибки численного интегрирования. Для расчета длинных траекторий был предложен алгоритм, основанный на идее квазилинвзриза. ции [140) и сохранении полной энергии системы вдоль траектории [49). Классический гамипьтониан системы и материальных точек, взаимодействующих друг с другом по закону, определяемому потенциальной энергией ((. в декартовой системе координат может быть записан в виде 1 з" Р. ог = — х — ~ + (.г (оы сг, -.. д ) .
зи (3.95) (3.97) 7В Обозначения выбраны здесь так, что три компоненты импульса и коор динаты каждой точки имеют индексы, которые пробегают последователь- ные значения. Например, первая материальная точка имеет компоненты импУльса Р>, Рт. Рз компоненты кооРдинаты — оы от, оз, и масса этой точки равна лт1 тт = тз. При численном решении будет использовать ся только декартова система координат.
Обычно ППЗ задается как функ. ция парных расстояний между точками. Эти расстояния однозначно опре деляются декартовыми координатами, а ориентация начальной конфигуре ции системы частиц относительно осей координат несущественна и вы. бирается произвольно, Уравнения движения Гамильтона для гамильтониана (3.95) имеют вид: ди Р- = г да, (3.96) Рг — Р (О! =Р .. С (О) =о ., )=1,...,3п.
,„' г ор г о' г Систему уравнений (3.96! можно записать в векторном виде: с(х — = 1(х), х (О) = хо, ггг где = [Рн Рт,".,Рзи.с1 ° Фг. -'Чзи) д(.г д(7 р, р „ (3.98) дд, д|! „гл, т Т о [РоыРог.-"Рози Чоычот. -.чоз„ вЂ” векторы-столбцы размерности (6Хл!. ЗИ р — — г+ сг, (ц,р). (3.99) 2 '=~ тг Вектор свободных параметров может быть найден при решении апгеб. раических уравнений, описывающих законы сохранения адднтивных инте градов движенигп 7т(р(д(г)). ц(д(г))) =(уе. йй(р(д(г) ).
а(и В) ) ) = м,, РВ(и(г)) =Р,, (3.100) 79 Используемые в расчетах ППЭ обладают непрерывными частными производными, что обеспечивает существование и единственность решения задачи Коши (3.97) для консервативной системы на любом отрезке времени. Рассмотрим свойства решений гамипьтоновых систем уравнений движения сложных молекул. Эти решения, как правило, осцилляционные, имеют разброс характеристических времен и подчиняются законам сохранения аддитивных интегралов движения, что предъявляет строгие требования к процедуре чиспенного расчета.
Известно (174), что поведение ошибки численного решения задачи Коши определяется спек~ром матрицы Якоби Л(х) = М/Ох. Если у матрицы 3 (х) действительная часть собственных значений положительна„ то с ростом времени растет и норма ошибки, т.е. решение системы неустой чиво. В случае отрицательной действительной чести собственных значений норма ошибки уменьшается и решение устойчиво. Прн наличии чисто мнимых собственных значений норма ошибки, возникающая при численном интегрировании, не убывает, что приводит к ее накоплению. Уравнения движения для консервативных систем имеют в основном мнимые собственные значения матрицы Якоби, что и является причиной осциппя. ционного характера решений. Зто обусловливает строгие требования к контроле точности численного решения.
Разброс характеристических времен делает задачу жесткой, что должно быть учтено при конструировании численного алгоритма. В работе (387) показано, что процедура численного решения классических уравнений движения должна обпадать свойствами консерватив ности. Несохранение аддитивных интегралов движения при чиспенном решении может приводить к неконтролируемому накоппенню ошибок. Предлагаемый алгоритм численного решения системы дифференциапь ных уравнений основан на методе локальной пинеаризацин (140). На каждом шаге интегрирования исходная ППЗ аппроксимируется квадрэ тичной формой, возниквощая при этом новая система дифференцнапь. ных уравнений является линейной н, спедоватепьно, допускает точное решение.
улучшая аппроксимацию, можно добиваться сходнмости нового решения к решению исходной задачи на всем интервале интегрированить Так как близкие поверхности определяют практически одинаковые модепи, то в смысле "траекторной нормы" решения должны сходиться.