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

Восстановление функции электронной плотности преобразованиями Фурье

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

Модель для компьютерного эксперимента:

На отрезке [0,30] (Å) расположены две молекулы. Атомы в молекуле связаны ковалентно и находятся на расстоянии 1 Å друг от друга. Молекулы расположены на расстоянии 5 Å. Всего 5 атомов. Электронные плотности (ЭП) атомов описываются гауссовой кривой. Максимум ЭП в центре атома приблизительно пропорционален числу электронов в атоме.
Гауссова функция определяется числами lambda,beta,gamma по формуле:
gauss = lambda*exp(-(beta^2)*(X-gamma)^2).
Функцию задали с помощью скрипта compile-func.py с соответствующими параметрами:
python compile-func.py -g 40,3,4+35,3.5,6+2,3,8+9,3,15+28,3,20
Получилась сумма пяти гауссовых функций с max в точках 4, 6, 8, 15, 20 и высотой 32, 14, 2, 12 и 16, соответственно.
На выходе имеем текстовый файл func.txt с парами значений X и Y и график функции:

Рис.1 График заданной функции ЭП

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

Коэффициенты Фурье рассчитали с помощью скрипта func2img.py На вход дается функция в файле func.txt. Выходной файл func_ft.txt содержит амплитуды и фазы гармоник с соответствующими номерами. Получилось 499 гармоник
Теперь проведем фильтрацию с помощью скрипта img-filter.py и найдем n_0, при котором восстановление отличное (т.е., при котором по графику восстановленной функции можно определить положение максимума всех гауссовых слагаемых функции ("атомов")).
Ниже приведены результаты восстановления по полным наборам для четырех n: 15, 25, 30, 34. n_0 определено как 34.

Рис.2 Графики восстановленной и исходной функций (слева; пунктирная линия соответствует восстановленной функции, сплошная - исходной) и только восстановленной (справа) функции при n=15.

Рис.3 n=25.

Рис.4 n=30.

Рис.5 n=34.

При n=15 и 25 восстановление можно назвать "средним": хорошо определяется положение четырех атомов. При n=30 восстановление "хорошее": можно было бы угадать положение и самого низкого пика, если бы было известно количество атомов. При n=34 восстановление можно назвать "отличным", определяются положения всех атомов.

Далее был добавлен шум к амплитудам и/или фазам при восстановлении по полному набору гармоник 0,...,34: 1) 10% по φ, 2) 10% по F, 3) 15% по φ и F, 4) 20% по φ и F. Результаты приведены ниже.

Рис.6 Графики восстановленной (по полному набору 0,...,n_0) функции без шума и с шумом по фазе 10% (слева) и только с шумом (справа).

Рис.7 Шум по амплитуде 10%.

Рис.8 Шум и по амплитуде, и по фазе 15%.

Рис.9 Шум и по амплитуде, и по фазе 35%.

Сравнив рисунки 6 и 7, можно заметить, что при добавлении шума 10% по фазе, колебания графика увеличились сильнее, чем при добавлении такого же шума по амплитуде. Однако атом водорода лучше угадывается на рисунке 6.
При шуме 10% восстановление можно считать "хорошим", при 15% - "средним", при 35% - "плохим".

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

Для получения неполного набора сначала удалим 1) гармонику 0, 2) гармоники 0 и 1. Для этого используем скрипт img-filter.pyс параметром -r 1-34 или -r 2-34. Результаты восстановления приведены ниже.

Рис.10 Пунктирная линия соответствует графику восстановленной по неполному набору гармоник 1-34 функции, сплошная линия - исходной функции.

Рис.11 Восстановление по неполному набору 2-34.

При исключении гармоники 0, график только смещается по вертикальной оси (гармоника 0 соответствует первому коэффициенту в ряду Фурье). При исключении двух первых гармоник функция искажается, но положения атомов по-прежнему четко определяются.

Теперь попробуем удалить 5-10% гармоник из середины набора: сначала удалим гармоники 17-19. Для этого используем img-filter.pyс параметром -r 0-16, 20-34. Результаты представлены ниже.

Рис.12 Восстановление по неполному набору, без гармоник 17-19.

Попробуем удалить побольше гармоник из середины, с 16 по 20.

Рис.13 Восстановление по неполному набору, без гармоник 16-20.

Получившиеся восстановления можно назвать "средними". Но чем больше гармоник удаляется из середины набора, тем сильнее максимумы от атомов похожи на шум и тем ниже качество восстановления.

Теперь добавим одну гармонику с номером, превышающим n_0 на 10 - гармонику 44 (img-filter.py -r 0-34,44).

Рис.14 Графики исходной и восстановленной функций при добавлении гармоники 44 к полному набору (слева) и без добавления (справа).

Значительных отличий от восстановления по полному набору 0,...,n_0 я не заметила, лишь слабо различимое увеличение шума.

Как определять разрешение для набора гармоник Фурье, по которым восстанавливается функция?

Разрешением полного набора гармоник называется период гармоники с наибольшим номером, а период гармоники равен расстоянию между соседними максимумами синусоиды.
В моем случае разрешение (d=T/n) составит 30/34=0,88 Å.

Полнота данных определяется как процент присутствующих в наборе гармоник с длиной волны, большей d, от максимально возможного при данном d числа гармоник. Для полного набора данных полнота равна 100%.

Для неполного набора нет строгого определения разрешения, кроме него необходимо также сообщать полноту данных. Так, например, для набора (0-15)+(21-34) разрешение равно 0,88 Å при полноте данных 85,7%. Так же можно рассуждать и для набора (2-34).
Однако для того случая, когда рассматривается набор (0-34)+44, такая логика не совсем верна. Получается, что разрешение равно 0,67 Å при полноте данных 80%. В таком случае логично было бы не рассматривать гармонику 44 и считать разрешение равным 0,88 Å при полноте 100%, и хуже от этого не станет.

Т.е., при определении разрешения стоит обращать внимание на структуру набора, и при малой полноте иногда можно выкинуть единичные гармоники, по номеру далеко стоящие от непрерывного или почти непрерывного ряда.

Таблица 1 Восстановление функции по коэффициентам ряда Фурье

Набор гармоник Разрешение (Å) Полнота данных (%) Шум амплитуды (% от F) Шум фазы (% от φ) Качество восстановления

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

0-15 2 100 0 0 среднее

0-25 1,2 100 0 0 среднее

0-30 1 100 0 0 хорошее

0-34 0,88 100 0 0 отличное

0-34 0,88 100 10 0 хорошее

0-34 0,88 100 0 10 хорошее

0-34 0,88 100 15 15 среднее

0-34 0,88 100 35 35 плохое

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

1-34 0,88 97 0 0 отличное

2-34 0,88 94 0 0 отличное

(0-16)+(20-34) 0,88 91 0 0 среднее

(0-15)+(21-34) 0,88 85,7 0 0 среднее

(0-34)+44 0,88 100 0 0 отличное

Набор гармоник	Разрешение (Å)	Полнота данных (%)	Шум амплитуды (% от F)	Шум фазы (% от φ)	Качество восстановления
Полный набор гармоник
0-15	2	100	0	0	среднее
0-25	1,2	100	0	0	среднее
0-30	1	100	0	0	хорошее
0-34	0,88	100	0	0	отличное
0-34	0,88	100	10	0	хорошее
0-34	0,88	100	0	10	хорошее
0-34	0,88	100	15	15	среднее
0-34	0,88	100	35	35	плохое
Неполный набор гармоник
1-34	0,88	97	0	0	отличное
2-34	0,88	94	0	0	отличное
(0-16)+(20-34)	0,88	91	0	0	среднее
(0-15)+(21-34)	0,88	85,7	0	0	среднее
(0-34)+44	0,88	100	0	0	отличное