Электронная плотность (одномерная модель), ее разложение в ряд Фурье и восстановление

Задание функции

Модель электронной плотности (ЭП) для компьютерного эксперимента

• Пусть имеются две молекулы на отрезке [0,30] ангстрем.
  
• Атомы в молекуле связаны ковалентно.
  
• Первая молекула состоит из двух атомов, между которыми расстояние в 1.25 ангстрем, а вторая из трех, между которвми расстояния равны 1.1 и 1.5
  
• Пусть молекулы расположены на расстоянии 4 ангстрем.
  
• Электронные плотности (ЭП) атомов описываются гауссовой кривой.
  
• Максимум ЭП в центре атома приблизительно пропорционален числу электронов в атоме.
  

Гауссова функция определяется числами lambda,beta,gamma по формуле:

gauss = lambda*exp(-(beta^2)*(X-gamma)^2)

.

Функция была создана с помощью скрипта compile-func.py суммой гауссовых функций, где первые позиции задают параметры высот колокола (высота зависит от числа электронов в атоме), вторые позции - ширина колоколов (порядка 1), третьи позиции - точки, в которых функции имееют максимумы.

Задание функции

По условию, молекулы состоят из разных атомов, что влечет за собой разницу в максимальных значениях их ЭП.
C помощью скрипта compile-func.py была задана функция со следующими параметрами:

• Были выбраны max в точках 5,6.25; 10.25, 11.35, 12.5.
• Высоты пиков - 40, 30, 2; 30, 20.
• Ширина колокола у всех атомов была взята порядка 1.
• Самый низкий (высота 2) пик моделирует атом водорода.

Команда:python compile-func.py -g 40,3,5+30,2.5,6.25+2,2,10.25+30,4,11.35+20,3,12.5

Был визуализирован суммарный график гауссовых функций, представленный на рисунке 1.

Рисунок 1. Суммарный график гауссовых функций.

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

С помощью скрипта func2fourier.py рассчитываются коэффициенты Фурье (на вход подается результат работы compile-func.py).
В результате был получен файл, содержащий 499 гармоник.
Команда:python func2fourier.py -i func.txt -o func_ft.txt

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

Для определения минимального числа гармоник, необходимых для отличного восстановления исходной функции ЭП, использовался скрипт fourier-filter.py.

Команда: python fourier-filter.py -i func_ft.txt -r 0-5 -o func_ft_filtered_0-5.txt
Для последующей визуализации восстановленной функции по отобранным гармоникам использовался скрипт fourier2func.py
В таблице 1 приведены результаты для разного количества выбранных полных наборов гармоник без шума.

Таблица 1. Результаты восстановления ЭП для разного количества выбранных полных наборов гармоник (0-5 и т.д.) 0-5 0-15 0-25 0-35 0-45 0-50 0-55 0-60

На 55 гармониках функция плотности атома водорода была восстановлена с хорошим качеством.
На 45-50 гармониках функция плотности атома водорода приобретает характерный сигнал, но для уточнения функции необходима дополнительная информация.
Видно, что для полного восстановления функции нужно примерно 60 гармоник.

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

Следующим шагом было добавление шума к амплитудам (параметр -F) и фазам (параметр -P) с помощью скрипта func2fourier.py. Отбирались первые 55 гармоник.
Это является лучшим приближением, когда удается получить лишь какой-то набор гармоник.


Результаты представлены в таблице 2.

Таблица 2. Результаты восстановления ЭП для разного количества выбранных полных наборов гармоник (0-5 и т.д.) с добавлением амплитудных и фазовых шумов. Амплитудный шум 20% 50% 70% Фазовый шум 20% 50% 70% Амплитудный шум 20%. Варьируется фазовый шум 20% 50% 70% Фазовый шум 20%. Варьируется амплитудный шум 20% 50% 70%

При добавлении шума самые значительные изменения вызвал фазовый шум.
При значениях шума 50-70% можно потерять информацию о средних по размеру атомах, взятых для симуляции.
При 70% амплитудного шума теряется информация о положении атома водорода у второй молекулы.
При комбинации самым "страшным" был вариант с комбинацией 20% амплитудного шума и 70% фазового шума. Мы потеряли почти все атомы.
Чтобы получить информацию об устойчивости восстановленной функции в зависимости от набра гармоник, было решено использовать следующие наборы гармоник:

• Одну-две начальные гармоники (номер 0, 0-1)
• Удалить 5-10% гармоник из середины набора (55*0.1 -> 23-28, 24-27)
• Добавить одну гармонику с номером, превышающим n_0 на 10 (65)

Для моделирования был использован набор гармоник из предыдущего эксперимента, который видоизменялся в соответствии с указанными правилами. Результаты представлены в таблице 3.

Таблица 3. Результаты восстановления ЭП для неполного набора гармоник -0 -(0-1) -(23-28) -(24-27) +65

Исходя из наблюдений можно сделать вывод, что разные гармоники вносят различный вклад.
Первые гармоники задают общее положение функции: при их удалении функция начинает опускаться и приобретать побочные синусоидальные колебания.
Гармоник из центра задают амплитуду высоты пиков и амплитуду функции. При удалении этих гармоник наблюдается общее снижение высот пиков и возникновение частых колебаний, приводящих к невозможности идентифицировать средние и малые атомы.
Интересно, что при добавлении гармоники №65 наблюдается ухудшение показателей функции по сравнению с набором 0-55. Думаю, это связано с переоцениванием веса этой гармоники. Когда мы делаем синтез Фурье, каждая последующая гармоника поправляет предыдущие и уточняет общий результат. Но если убрать промежуточные поправки, то гармоника №65 может привнести не вполне адекватную информацию.
Я попробовал добавить больше гармоник с большими номерами (65-69), но это не помогло. Видимо, предыдущие гармоники более важны.

Рисунок 2. Синтез Фурье для гармоник 0-55 + 65-69 .

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

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

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

1. PDB: 3R1R
2. RCSB: 3R1R