• На главную
• Семестры
• Проекты
• Заметки
• Обо мне
• Ссылки

Восстановление функции электронной плотности

Цель и задачи

Цель данной работы заключается в исследовании качества восстановления функции ЭП от одной переменной в зависимости от того, какие и сколько гармоник ряда Фурье используются для ее восстановления.

Задание включает следующие этапы:

1. Создание модельной функции ЭП в одномерной элементарной ячейке
2. Расчет параметров сигнала, моделирующих экспериментальные данные: амплитуды и фазы
3. Восстановление функции ЭП по модельным ("экспериментальным") данным
4. Оценка качества восстановления функции ЭП

Создание модельной функции ЭП в одномерной элементарной ячейке

Была создана линейная атомная модель для компьютерного эксперимента: 2 линейных молекулы O-C-O и Cl-Hg-Cl. Для этой модели была смоделирована функция, описывающая одномерную электронную плотность (ЭП) с помощью скрипта compile-func.py:

ЭП атомов описываются гауссовской кривой, максимум в центре атома приблизительно пропорционален числу электронов в атоме. Функция задается на интервале [0Å, 30Å]. Гауссова функция определяется числами lambda, beta, gamma по формуле:

lambda - высота гауссовской кривой, соответствующая максимуму электронной плотности в центре атома. Пропорционально числу электронов в атоме:

• O - 32
• C - 12
• Hg - 201
• Cl - 71

beta - ширина колокола, которая примерно соответствует диаметру атома:

• O - 1.46
• C - 1.82
• Hg - 2.98
• Cl - 2.04

gamma - координата максимума гауссовой кривой. Расстояние между молекулами - 5.5 Å, длина всех связей - 1.5 Å. График электронной плотности модели атомов представлен на рисунке 1. Поточенчные координаты приведены в файле func.txt.

Рис.1 График электронной плотности 2 линейных молекул O-C-O и Cl-Hg-Cl.
По оси абсцисс отложены кординаты атомов в Å, а по оси ординат отложена амплитуда ЭП в условных единицах.

Расчет амплитуд и фаз сигналов, моделирующих экспериментальные данные

Рассчет коэффициентов Фурье

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

Коэффициенты Фурье рассчитываются с помощью скрипта func2fourier.py следующей командой:

Скрипт по ранее полученному файлу func.txt выдает 499 гармоник на отрезке [0 Å, 30 Å]. Выходной файл скрипта – файл func_ft.txt. Выходной файл имеет формат: <номер гармоники> <амплитуда> <фаза>.

Отбор гармоник

Набор гармоник ряда Фурье называется полным, если известны все гармоники с номерами 0, 1, 2, ..., n. Разрешением полного набора гармоник называется период гармоники с номером n, то есть с наибольшим номером. Период гармоники равен расстоянию между соседними максимумами синусоиды. Его также называют длиной волны этой гармоники. Зависимость длины волны гармоники от номера:

• 0 – гармоника есть константа, длины волны нет
• 1 – длина волны равна длине отрезка T (в нашем случае T = 30 Å)
• 2 – — " — T/2
• ...............
• n – — " — T/n

Отбор гармоник для получения неполного набора можно сделать одним из способов:

• Скопировать файл с коэффициентами Фурье и убрать лишние строки.
• Пометить лишние строки знаком "#" в первой позиции строки (такие строки при чтении игнорируются).
• Запустить скрипт fourier-filter.py, который отфильтрует нужные строки.

Мы будем пользоваться последним способом.

Восстановление функции ЭП по амплитудам и фазам части сигналов

Полные наборы гармоник

По полному набору гармоник был построен график восстановленной функции электронной плотности модели с помощью скрипта fourier2func.py командой:

Сравнение графика восстановленной функции и графика первоначальной функции можно увидеть на рис.2. Очевидно, что при количестве гармоник, равном 499, получается отличное восстановление, т.е. по графику восстановленной функции можно определить положение максимумов всех гауссовых слагаемых. При таком количестве гармоник разрешение составляет ... Å (разрешением полного набора гармоник называется период гармоники с номером n, тo есть с наибольшим номером. Период гармоники равен расстоянию между соседними максимумами синусоиды). Это довольно высокое значение.

Рис.2 График электронной плотности 2 линейных молекул O-C-O и Cl-Hg-Cl
По оси абсцисс отложены кординаты атомов в Å, а по оси ординат отложена амплитуда ЭП в условных единицах.

A. График первоначальной функции	B. График восстановленной функции

Найдем минимальное число гармоник, при котором восстановление функции сохранится отличным. Для этого рассмотри n-первых гармоник, где n = 2, 21, 26, 31, 41.

Как сравнить восстановленную функцию с исходной?

• Отличное восстановление – по графику восстановленной функции можно определить положение максимума всех гауссовых слагаемых функции ("атомов");
• Хорошее восстановление – можно угадать положение всех максимумов, зная число слагаемых ("атомов"), хотя на восстановленной функции максимумы от атомов не отличимы от шума;
• Среднее восстановление – положение каких-то атомов определить по восстановленной функции нельзя, других - можно;
• Плохое восстановление – положение атомов определить не представляется возможным; можно только предсказать примерный размер "молекулы".

Рис.3 Графики восстановленных функций
Графики восстановленных функций, представленны красной линией, для A. - 0-1; Б. - 0-5; В. - 0-10; Г. - 0-25; Д. - 0-20 и Е. - 0-30 наборов гармоник. Исходная функция – представлена в виде сплошной линии.

Исходя из графиков на рис.2, видно, что 31 - это минимальное число гармоник (n₀), при котором функция сохраняет отличное восстановление.

Добавим шум к амплитудам и фазам при восстановлении по полному набору гармоник. Добавление гауссового шума к амплитудам (параметр -F <число>) и фазам (-P <число>) искажает все вычисленные амплитуды и фазы. Пример: -F 20 (шум 20%) приводит к тому, что к каждой амплитуде прибавляется случайное число, распределенное нормально с параметрами: среднее=0, среднее квадратичное отклонение σ=0.2*F. Аналогично действует параметр -P <число>. Добавим 10% шума только к амплитудам, только к фазам, по 5 % к амлитудам и фазам:

Рис.4 Графики восстановленных функций
Графики изначальной функции электронной плотности (черная линия) и восстановленной (красная линия) для 31 гармоник A. - с шумом 10% к амплитуде; Б. - 10% к фазе; В - по 5% к амплитуде и фазе.

Неполные наборы гармоник

Сделаем неполные наборы гармоник для набора n₀: удалим несколько первых гармоник, удалим серединные гармоники, добавим гармонику с номером n₀+10.

При удалении нулевой гармоники функция заметно смещается вниз по оси ординат. Это можно объяснить тем, что эта гармоника есть константа, длины волны у нее нет. Если же удалить первые три, то вычитается синусоида с периодом T/3. Видно, что в этом случае мы наблюдаем "провал" амплитуды для средних значений координат.

Рис.4 Графики восстановленных функций
Графики восстановленных функций, представленны красной линией, для A. - Удаление первой гармоники; Б. - Удаление первых трех гармоник; В. - Удаление 5-10% гармоник из середины набора; Г. - Добавление гармоники с номером n₀+10. Исходная функция – представлена в черной линией.

Вследствие удаления 10% гармоник в середине набора качество восстановления сильно ухудшилось, максимумы от атомов стали почти не отличимы от шума. Добавление гармоники с номером 40 к набору гармоник n₀ немного "дестабилизировали" шум, по сравнению с просто набором гармоник n₀. То есть шум стал как бы неравномерным, а где-то меньше, где-то больше. Из всего выше сказанно можно сделать вывод, что чем больше номер гармоники, тем меньший вклад она вносит в общий вид восстановленной функции ЭП. Это логично, так как с увеличением номерам период гармоники уменьшается.

Cравнение восстановленной функцию с исходной

Для неполного набора данных нет строгого определения разрешения. Кроме разрешения d необходимо сообщить полноту данных — процент гармоник с длиной волны большей d от максимально возможного, присутствующих в наборе. Для полного набора данных (разрешение d=T/n) полнота равна 100%. Рассмотрим алгоритм поиска разрешения для неполного набора гармоник (0-30, 40). Если применить разрешение для полного набора, то получится, что разрешение = 0.75, а полнота данных всего 80% (=32/40). Поэтому разумнее ввести порог на полноту данных 90% и пересчитать наше разрешение.

© Nuzhdina Ekaterina, 2015