import numpy as np
import pandas as pd
import prody as prd
import matplotlib

ls 

02.txt              30.txt              Prak4.1.pse         func2fourier.py
04.txt              40.txt              Site_prak4.pages*   lections.txt
05.txt              50.txt              Task1.png           mol-replacement.py
08.txt              60.txt              Task3((((.ipynb     out.txt
10.txt              70.txt              __pycache__/        outF25.txt
12.txt              7bk2.pdb.gz         combine-show.py     outF50.txt
14.txt              F25_60.txt          combine.py          outF75.txt
15.txt              F50_60.txt          compile-func.py     outP25.txt
16.txt              F75_60.txt          fourier-filter.py   outP50.txt
18.txt              P25_60.txt          fourier.py          outP75.txt
20.txt              P50_60.txt          fourier2func.py     two_func.txt
22.txt              P75_60.txt          func.txt

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

Рассмотрим одномерную вселенную, в которой атомы находятся на одном отрезке длиной 30 ангстрем. Сгенерируем функцию электронной плотности на интервале [0,30] с помощью скрипта compile-func.py. Функция будет иметь вид суммы нескольких гауссовых кривых с центром в разных точках. При этом, центры атомов, находящиеся на расстоянии 1-1.5 ангстрем считаются ковалентно связанными, а расстояние 3-5 ангстрем свойственно для отдельных молекул. Гауссова функция определяется числами lambda, beta, gamma по формуле: gauss = lambdaexp(-(beta^2)(X-gamma)^2) Здесь a^2=a*a, где Lambda задает высоту гауссиана, через которую мы имитируем число электронов у атома, Beta задает ширину гауссиана, значение 3 соответствует релевантной нам ширине колокола плотности около 1 ангстрем, что похоже на реальную ситуацию, Gamma задает положение центра атома. [https://kodomo.fbb.msu.ru/wiki/2018/7/task4]

Модель электронной плотности одномерной системы атомов из двух молекул, состоящих из 2 и 3 атомов соответственно:

run compile-func.py -g 30,3,3.2+20,3,4.4+40,3,5.4+15,3.5,7.8+2,3,9

File func.txt created

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

Теперь смоделируем получение экспериментальных даннных (амплитуды и фазы), которые рассчитываются однозначно по входной функции ЭП с помощью скрипта func2fourier.py (он вычисляет коэффициенты Фурье).

run func2fourier.py -i func.txt -o out.txt 

Done
Done

#исправляю print"..Done" на print("Done"), убираю многочисленные "\"

Выходной файл out.txt имеет формат <номер гармоники> <амплитуда> <фаза>.

Описанние дальнейших действий и цели практикума

В эксперименте амплитуды опропределяются не для всех сигналов. Кроме того, интенсивности сигналов (амплитуды) и фазы определяются с ошибкой. Поэтому далее я буду "портить" полученные идеальные экспериментальные данные следующими способами: 1) Буду использовать небольшой интервал полного набора гармоник (с помощью скрипта fourier-filter.py, который отфильрует нужные строки)

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

3) Буду использовать неполный набора гармоник (с помощью скрипта fourier-filter.py, который отфильтрует нужные строки) Результаты буду заносить в ТАБЛИЦУ

В процессе действи 1-3 буду анализировать, как именно неполнота данных из эксперимента влияет на качество восстановления электронной плотности, сравнивая ее с исходной (идеальной) . Для классификации восстановления по качеству были предложены следующие критерии:

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

Сам Анализ

Использую разные интервалы полного набора гармоник

Рассмотрим несколько полных наборов гармоник без добавления шума. Буду использовать следующие интервалы: 0-2, 0-10, 0-15, 0-20, 0-30, 0-40, 0-50, 0-60

run fourier-filter.py -i out.txt -r 0-2 -o 02.txt

..Done

Восстановление функции по отобранным гармоникам (вспомните, что такое гармоника) выполняется скриптом fourier2func.py

run fourier2func.py -f func.txt -i 02.txt 

File two_func.txt with function and its recovering is created

Качество восстановления плохое. Мы даже не можем определить количество молекул в структуре, только их примерное положение

run fourier-filter.py -i out.txt -r 0-10 -o 10.txt

..Done

run fourier2func.py -f func.txt -i 10.txt 

File two_func.txt with function and its recovering is created

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

run fourier-filter.py -i out.txt -r 0-15 -o 15.txt

..Done

run fourier2func.py -f func.txt -i 15.txt

File two_func.txt with function and its recovering is created

Качество электронной плотности резко возрастает до среднего. Можем определить положение трех атомов из пяти

run fourier-filter.py -i out.txt -r 0-20 -o 20.txt

..Done

run fourier2func.py -f func.txt -i 20.txt

File two_func.txt with function and its recovering is created

Качество для первой молекулы - среднее, а для второй уже хорошее. Несмотря на то, что на восстановленной функции максимум от водорода не отличимы от шума, можно угадать его положение, зная атомов.

run fourier-filter.py -i out.txt -r 0-30 -o 30.txt 

..Done

run fourier2func.py -f func.txt -i 30.txt 

File two_func.txt with function and its recovering is created

run fourier2func.py -f func.txt -i 30.txt -s

File two_func.txt with function and its recovering is created

Качество молекулы, состоящей из трех атомов, теперь можно назвать отличным! Восстановленная электронной плотность совпадает с исходной. Однако качество восстановления для второй молекулы я бы по прежнему назвала хорошим (видно, если воспользоваться параметром -s, центр водорода смещен)

run fourier-filter.py -i out.txt -r 0-40 -o 40.txt 

..Done

run fourier2func.py -f func.txt -i 40.txt 

File two_func.txt with function and its recovering is created

run fourier-filter.py -i out.txt -r 0-50 -o 50.txt

..Done

run fourier2func.py -f func.txt -i 50.txt -s

File two_func.txt with function and its recovering is created

Качество для второй молекулы пока можно назвать только хорошим, положение водорода нельзя определить точно так как "колокол" Гауссовой функции несимметричен

run fourier-filter.py -i out.txt -r 0-60 -o 60.txt

..Done

run fourier2func.py -f func.txt -i 60.txt

File two_func.txt with function and its recovering is created

А вот теперь качесвто функции электронной плотности отличное для обеих молекул

Качество восстановления функции электронной плотности увеличивается с увеличением длины интервала гармоник, так как мы с большей точностью способы апроксимировать исходную функцию.

Добавлю шум амплитуды - 25%, 50%, 75%; буду использовать при этом интервал полного набора гармоник, при котором наблюдается отличное восстановление ЭП (0-60)

run func2fourier.py -i func.txt -o outF25.txt -F 25

Done
Done

run fourier-filter.py -i outF25.txt -r 0-60 -o F25_60.txt

..Done

run fourier2func.py -f func.txt -i F25_60.txt 

File two_func.txt with function and its recovering is created

При 25% уровне зашумления качество певой молекулы всё ещё можно назвать отличным, но качество второй молекулы портится до хорошего. Пик для водорода можно принять за шум, если не знать, что он есть в составе.

run func2fourier.py -i func.txt -o outF50.txt -F 50

Done
Done

run fourier-filter.py -i outF50.txt -r 0-60 -o F50_60.txt

..Done

run fourier2func.py -ffunc.txt -i F50_60.txt

File two_func.txt with function and its recovering is created

Мы все еще можем определить точное положение атомов в первой структуре - для нее качестов отличное. А качество для второй молекулы снижается до плохого - атом водорода неразличим.

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

run func2fourier.py -i func.txt -o outF75.txt -F 75

Done
Done

run fourier-filter.py -i outF75.txt -r 0-60 -o F75_60.txt

..Done

run fourier2func.py -f func.txt -i F75_60.txt 

File two_func.txt with function and its recovering is created

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

Добавлю шум фаз - 25%, 50%, 75%; буду использовать при этом интервал полного набора гармоник, при котором наблюдается отличное восстановление ЭП (0-60)

run func2fourier.py -i func.txt -o outP25.txt -P 25

Done
Done

run fourier-filter.py -i outP25.txt -r 0-60 -o P25_60.txt

..Done

run fourier2func.py -f func.txt -i P25_60.txt

File two_func.txt with function and its recovering is created

Здесь шум фазы сосатвляет 25%. Качество первой молекулы отличное, но качество второй молекулы портится до хорошего, как и в случае добавления 25% шума к амплитудам.

run func2fourier.py -i func.txt -o outP50.txt -P 50

Done
Done

run fourier-filter.py -i outP50.txt -r 0-60 -o P50_60.txt

..Done

run fourier2func.py -f func.txt -i P50_60.txt

File two_func.txt with function and its recovering is created

Точное положение атомов в первой структуре без знания ее состава определить не получится - для нее качество хорошее, близкое к среднему. А качество для второй молекулы снижается до плохого - атом водорода неразличим. Кроме того, шум настольк велик, что можно предположить, что вторая молекула находится сильно правее. Этот результат хуже, чем тот, который был получен при шуме амплитуды 50%, что соответствует теоритическим знаням о бОльшем вкладе фазы в определения элетронной плотности.

run func2fourier.py -i func.txt -o outP75.txt -P 75

Done
Done

run fourier-filter.py -i outP75.txt -r 0-60 -o P75_60.txt

..Done

run fourier2func.py -f func.txt -i P75_60.txt

File two_func.txt with function and its recovering is created

Зашумление слишком сильное, чтобы сказать про структуры хоть что-нибудь. С ней совсем нельзя работать. Зашумление фаз очень сильно повлияло на итоговые данные.

Скомбинирую шумы (30% для фаз, 10% для амплитуды и наоборот), буду использовать при этом интервал полного набора гармоник, при котором наблюдается отличное восстановление ЭП (0-60)

run func2fourier.py -i func.txt -o combine1.txt -P 30 -F 10

Done
Done

run fourier-filter.py -i combine1.txt -r 0-60 -o combine1_060.txt

..Done

run fourier2func.py -f func.txt -i combine1_060.txt

File two_func.txt with function and its recovering is created

run func2fourier.py -i func.txt -o combine2.txt  -P 10 -F 30

Done
Done

run fourier-filter.py -i combine2.txt -r 0-60 -o combine2_060.txt

..Done

run fourier2func.py -f func.txt -i combine2_060.txt 

File two_func.txt with function and its recovering is created

Комбинация шумов, несмотря на небольшие значения (30 и 10%), очень сильно портит качество функции электронной плотности. Причем вновь при бОльшем шуме фаз, качество хуже и пики от шума выше.

Использую разные интервалы (2-60, 4-60, 8-60, 16-60) неполного набора гармоник

run fourier-filter.py -i out.txt -r 2-60 -o 260.txt

..Done

run fourier2func.py -f func.txt -i 260.txt

File two_func.txt with function and its recovering is created

run fourier-filter.py -i out.txt -r 10-60 -o 1060.txt

..Done

run fourier2func.py -f func.txt -i 1060.txt 

File two_func.txt with function and its recovering is created

run fourier-filter.py -i out.txt -r 20-60 -o 2060.txt

..Done

run fourier2func.py -f func.txt -i 2060.txt 

File two_func.txt with function and its recovering is created

При снижении полноты данных постепенно снижается и качество функции электронной плотности. Однако я ожидала, что качество даже для неполного набора гармоник 2-60 будет совсем плохим так как вклад первых гармоник в аппроксимацию функции электронной плотности наибольший.