Восстановление структуры по "экспериментальным данным PCA" в одномерной модели

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

Для компьютерного эксперимента была использована модель со следующими характеристиками:

• на отрезке [0,30] Å расположены две молекулы (молекула А из 2-х атомов и молекула В из 3-х атомов);
• атомы в молекуле связаны ковалентно и находятся на расстоянии 1-1.5 Å друг от друга;
• молекулы расположены на расстоянии 4.3 Å (водородная связь или гидрофобное взаимодействие между ними);
• электронные плотности (ЭП) атомов описываются гауссовой кривой, максимум ЭП в центре атома приблизительно пропорционален числу электронов в атоме.

Гауссова функция определяется числами λ, β, γ по формуле: gauss = λ*exp(-(β^2)*(X-γ)^2).

С помощью этой функции для компьютерной модели был получен 1D график электронной плотности. Для этого был использован скрипт compile-func.py, который был запущен с параметрами:

python compile-func.py -g 15,3,10.5+10,3.5,11.7+20,3.7,16+30,3.5,17.3+25,3,18.7

Каждая из троек чисел под параметром -g соответствует параметрам λ, β, γ в выше приведенной формуле, а число таких троек соответствует числу атомов в модели. При этом λ — высота пика (отражает число электронов в атоме), β — ширина, а γ — координата максимума пика.

В результате был получен файл func.txt с координатами точек и изображение смоделированной функции, представленное на Рис. 1.

Рис. 1. Изображение смоделированной функции плотности

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

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

Для разложение гауссовской функции в ряд Фурье был использован скрипт func2fourier.py, которому на вход подается полученный выше файл func.txt. В результате был получен файл func_ft.txt, содержащий 499 гармоник, для каждой из которых были восстановлены амплитуда (F) и фаза (phi).

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

Для отбора минимального числа гармоник, необходимых для хорошего восстановления исходной функции ЭП, были использованы скрипты fourier-filter.py и fourier2func.py. Первый скрипт использовался для отбора гармоник от 0 до n-й (10 < n < 50), а второй — для восстановления функции по отобранным гармоникам. Для каждого из 9 наборов гармоник была восстановлена ЭП с использованием команд (результаты на изображениях ниже):

python fourier-filter.py -r 0-n
python fourier2func.py -f func.txt -i func_ft_filtered.txt

Из изображений видно, что уже при n = 30 количество атомов и их положение определяются достаточно однозначно, но все же числом гармоник, нужным для отличного восстановления, было решено считать n = 40, т.к. тогда возможно точно определить максимумы каждого из пиков, а шумы практически исчезают.

Однако, все это верно только для идеального эксперимента, в котором отсутствуют ошибки в измерении интенсивности сигналов (амплитуд) и фаз. Поэтому далее к амплитудам (параметр -F <число>) и фазам (-P <число>) был искусственно введен гауссовский шум, который искажает все вычисленные амплитуды и фазы. Для этого были использованы команды, указанные ниже, а результаты представлены на изображениях ниже.

python func2fourier.py -i func.txt -F % -P %
python fourier-filter.py -r 0-40
python fourier2func.py -f func.txt -i func_ft_filtered.txt

При рассмотрении изображений видно, что при шуме по амплитудам 20% и по фазе 30% уже сложно определить положение атомов, а при шуме по амплитудам 30% — попросту невозможно. Аналогичная картинка видна при введении шума по обоим параметрам сразу — уже на уровнях шума 20% появляются ложные пики, сопоставимые по размерам с наименьшими истинными.

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

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

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

Сводные данные по качеству восстановления

Следует отметить "критерии качества" восстановления, используемые в данном задании:

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

С их использованием была составлена сводная таблица с информацией по качеству восстановления функции ЭП из различных рассмотренных наборов гармоник (Таблица 1).