Учебный сайт Макаровой Надежды
Восстановление функции электронной плотности преобразованиями Фурье
Задание функции. Модель.
В качестве примера для восстановления электронной плотности была использована следующая модель:
на отрезке [0,30] (ангстремы) расположены две молекулы (H-C-N и С-S-O). Атомы в молекуле связаны ковалентно и находятся на расстоянии 1- 1.5 анстрем друг от друга.
Молекулы расположены на расстоянии 3-5 ангстрем. Всего 6 атомов (3+3).
Электронные плотности (ЭП) атомов описываются гауссовой кривой. Максимум ЭП в центре атома приблизительно пропорционален числу электронов в атоме.две молекулы.
Функция электронной плотности атомов на этом отрезке была задана с помощью скрипта compile-func.py:
python compile-func.py -g 2,3,3+12,3,4.1+14,3,5.2+12,3,9.5+32,3,11+16,3,12
Первая цифра - высота пика (H-2, C-12, N-14, C-12, S-32, O-16)
Вторая цифра - ширина пика
Третья цифра - координата максимума
Выходной файл скрипта compile-func.py - func.txt, содержит для каждой точки на отрезке [0,30] значение функции в точке. Скрипт также позволяет получить изображение заданной совокупности гауссовых функций:
Полные наборы гармоник
Коэффициенты разложения функции в ряд Фурье были получены с помощью скрипта func2fourier.py:
python func2fourier.py -i func.txt
Выходной файл скрипта func2fourier.py – func_ft.txt. В нём приведены амплитуды и фазы гармоник, а также их номера.
Восстановление функции электронной плотности по набору гармоник
В первую очередь требовалось найти для 499 гармоник, полученных в предыдущем пункте, n0 - число гармоник, при котором восстановление отличное. Проверялись n0 =10, 20, 30, 40, 50. Ниже представлены графики, показывающие восстановленные функции электронной плотности. n0 было определенно равным 40.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Добавление шума
Реальные данные всегда содержат шум, поэтому скриптом func2fourier.py были созданы разложения функции ЭП, содержащие разный процент шума по амплитуде и/или фазе и включающие в себя полный набор из 40 гармоник. Процент шума задали как 10%, 20%, 30%.
| Без исходной | Без исходной | ||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Как видно на графиках при значениях 30% шума как по амлитуде, так и по фазе только опредление атомов водорода было затруднено. Однако шум по фазе значительно ухудшал восстановление, чем шум по амплитуде.
|
|
|
|
|
|
|
|
Внесение же шума (30%) и в фазу, и в амплитуду не позволяло определить "атомы углерода и азота".
Восстановление функции ЭП по неполному набору гармоник
Удаление начальных гармоник
Для получения неполного набора гармоник были удалены одна и две начальные гармоники (гармоника с номером 0; гармоники с номерами 0 и 1):
# Удаление гармоники c номером 0 python fourier-filter.py -r 1-40 -i func_ft.txt -o ft_1-40.txt python fourier2func.py -i ft_1-40.txt -o two_func_1-40.txt # Удаление гармоник с номером 0 и 1 python fourier-filter.py -r 2-40 -i func_ft.txt -o ft_2-40.txt python fourier2func.py -i ft_2-40.txt -o two_func_2-40.txt
Были получены следущие графики электронной плотности:
|
|
|
|
В эксперименте первые гармоники попадают в заглушку детектора и не могут быть определены. Однако их отсутвие в данном случае не сильно мешает определить максимумы, даже у "атомов водорода".
Удаление гармоник из середины набора
Ниже представлены графики, когда из полного набора (1-40) удалялись гармоники в первом случае (19-21) - 7.5%, во втором случае (16,17,21-25) - 17.5%
Уже на 7% трудно отличить пик "атома водорода" от шума. А при 17.5% и пик от "атома углерода". Добавление же гармоники с номером 50 никак не помогло улучшить ситуацию.
|
|
|
|
|
|
Разрешение
Разрешение для полного набора гармоник соответствует наименьшему из разрешений составляющих его гармоник (разрешение одной гармоники - период)
Вычисляется как длина рассматриваемого отрезка, деленная на номер последней гармоники (так как они упорядочены, и последняя имеет наименьший период):
d_0 = T/n . В рассматриваемом случае Т=30.
Полнота данных вычисляется как процент гармоник в наборе, период (= разрешение) которых не меньше разрешения, от числа гармоник в полном наборе для этого разрешения.
Для полного набора полнота данных равна 100%.
Для неполного набора гармоник определение разрешения отличается введением параметра полноты данных.
Если задаться полнотой, например, в 95%, то разрешение (d) неполного набора гармоник — это разрешение (d) такой гармоники,
что 95% гармоник из набора имеют разрешение больше d.
Выводы:
