Open: Обработка конформационных ансамблей в Python

New title: Обработка конформационных ансамблей в Python

В чем заключается эта обработка и для чего её делать?

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

• Отсеять конформации с чрезмерно высокой энергией;

• Удалить дубликаты;

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

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

• Объединить конформационные ансамбли из нескольких расчетов, полученные, например, как в пункте выше, и отсортировать структуры по их энергии.

Этот мануал посвящен моей библиотеке pyxyz, которая позволяет решать эти задачи несколькими строчками в Python. Также этот мануал может быть полезен при анализе XYZ-файлов, полученных в результате расчета молекулярной динамики.

Установка библиотеки

Для работы под Linux требуется Python версии 3.7 или выше. Для работы под Windows требуется Python версии 3.8 или выше.

Для установки на Linux нужно выполнить в терминале следующую команду:

pip install -i https://test.pypi.org/simple/ "pyxyz<1.0"

Для установки на Windows нужно выполнить в командной строке эту строчку:

pip install -i https://test.pypi.org/simple/ pyxyz>=1.0

Обратите внимание, что в некоторых случаях может понадобиться заменить pip на pip3 (это зависит от того, как у вас установлен Python).

После вызова pip install проверьте корректность установки: импортируется ли библиотека и выполняются ли тесты:

import pyxyz
import pyxyz.test
pyxyz.test.run_tests()

Пример обработки набора конформаций

Архив с примерами -- Pyxyz_manual.zip (File:Pyxyz manual.zip)

В качестве примера рассмотрим молекулу из файла diol_single.xyz и её конфнабор из CREST - diol_crest.xyz (см. архив). Рассмотрим несколько учебных задач, которые могли бы возникнуть на практике.

Слияние конфнаборов и сортировка

Допустим, мы знаем, что энергия структуры в файле diol_single.xyz равна -24.961, и мы добавить эту структуру к набору в файле diol_crest.xyz. Для начала, загружаем оба xyz-файла в один объект Confpool:

ensemble = Confpool()
ensemble.include_from_file("diol_crest.xyz", energy=lambda x: float(x.strip()))
ensemble.include_from_file("diol_single.xyz", energy=lambda x: -24.961)

Здесь мы передаем в include_from_file функцию energy, которая из описания структуры в xyz-файла вычленяет её энергию в виде дробного числа. Далее можно отсортировать полученный набор и сохранить его в файл:

ensemble.sort()
ensemble.save("result.xyz")

Анализ геометрии

Допустим, мы хотим узнать, сколько структур в наборе diol_crest.xyz имеют торсионный угол (5,6,7,19) близкий к 180. Для начала, готовим объект Confpool:

ensemble = Confpool()
ensemble.include_from_file("diol_crest.xyz", energy=lambda x: float(x.strip()))

И затем пользуемся функцией dihedral_count, чтобы узнать количество структур:

d = [5, 6, 7, 19]
n_conf = ensemble.dihedral_count(d[0], d[1], d[2], d[3], lambda x: abs(180.0 - x) > 10.0)
# OR n_conf = ensemble.dihedral_count(*d, lambda x: abs(180.0 - x) > 10.0)
print("Number of structures = {}".format(n_conf))

Требование на двугранный угол выражается функцией lambda x: abs(180.0 - x) > 10.0 - то есть отклонение от 180 должно быть не более 10 градусов. Обратите внимание, что распаковка листа на аргуметы делается короче, если пользоваться звёздочкой *.

Можно не только узнать количество структур, но и отфильтровать интересующие. Например, так выглядит код для фильтрования структур с контактами между парой атомов (1, 19) и сохранения результата в файл:

ensemble.distance_filter(1, 19, lambda x: x < 2.0)
ensemble.save("result2.xyz")

Передача данных в Python

Все данные, которые содержатся в объекте Confpool, могут быть

best_conf = ensemble.get_structure(0)
atom_symbols = ensemble.get_atom_symbols()
for i, symbol in enumerate(atom_symbols):
    print("%2s %14.6f %14.6f %14.6f" % (symbol, *best_conf['xyz'][i]))

Функционал библиотеки

Описание методов класса Confpool приведено по следующей схеме:

Имя функции(название аргумента: тип аргумента, ...) -> тип возвращаемого значения

Добавление структур в набор

include_from_file(filename: str, energy=None or function) -> None, energy(descr: str) -> float
Загрузка структур из XYZ-файле filename и (опционально) парсинг энергии из описания каждой структуры, с помощью функции energy.
filename - имя xyz-файла для импорта;
energy - keyword, который можно игнорировать или передавать функцию для выделения энергии из строчки описания в xyz-файле.
 
Пример, когда в описании к структуре энергия дается, как число без дополнительных подписей:
e.include_from_file("file.xyz", energy=lambda x: float(x.strip()))

Изменение данных

sort() -> None
Сортировка ансамбля по возрастанию энергии.
Пример: e.sort()

update_description(f: function) -> None, f(e: float, d: str) -> str
Используя функцию f сгенерировать новые описания к структурам из их энергии (float, первый аргумент f) и старого описания (str, второй аргумент f).
Пример: e.update_description(lambda ener, d: "E = {} a.u. (old descr = {})".format(ener, d))

Фильтры:

energy_filter(e_cutoff: float, etype='kcal/mol' or 'hartree'/'a.u.') -> None
Отбросить все конформеры, энергия которых выше энергии минимального на e_cutoff. В случае, когда etype=kcal/mol, разницы по энергии переводятся из a.u. в kcal/mol, после чего производится сравнение с e_cutoff. Во всех случаях предполагается, что энергия дана в a.u. Дефолтный etype=kcal/mol.
Пример: e.energy_filter(5.0, etype='kcal/mol')

distance_filter(idxA: int, idxB: int, condition: function) -> None, condition(dist: float) -> bool
Отбросить все конформеры, расстояние между атомами idxA и idxB которых не удовлетворяет условию condition. Функция condition получает на вход расстояние между атомами и, если она возвращает False, то структура отбрасывается.
Пример, в котором отбросятся структуры с расстоянием между атомами 1 и 2 больше, чем 1.7: e.distance_filter(1, 2, lambda x: x < 1.7)

valence_filter(idxA: int, idxB: int, idxC: int, condition: function) -> None, condition(angle: float) -> bool
Отбросить все конформеры, валентный угол между атомами idxA, idxB и idxC которых не удовлетворяет условию condition. Функция condition получает на вход угол между атомами в градусах и, если она возвращает False, то структура отбрасывается.
Пример, в котором отбросятся структуры в углами между атомами 1, 2 и 3 близкими к развернутому: e.valence_filter(1, 2, 3, lambda a: 180.0 - a < 10.0)

dihedral_filter(idxA: int, idxB: int, idxC: int, idxD: int, condition: function) -> None, condition(angle: float) -> bool
Отбросить все конформеры, торсионный угол между атомами idxA, idxB, idxC и idxD которых не удовлетворяет условию condition. Функция condition получает на вход угол между атомами в градусах (между -180 и 180) и, если она возвращает False, то структура отбрасывается.
Пример, в котором отбросятся структуры в торсионными углами между атомами 1, 2, 3 и 4 близкими к нулю: e.dihedral_filter(1, 2, 3, 4, lambda a: abs(a) < 10.0)

Получение данных

save(filename: str) -> None
Выписать данный конформационный набор в файл с название filename в формате XYZ. Струтуры будут выписаны в определенном порядке, например, в котором они находились в ранее импортированном файле или порядке после сортировки по энергии. В полученном XYZ файле будут включены описания к каждой структуре.
Пример: e.save("result.xyz")

size() -> int
Возвращает количество структур в данном конформационном наборе.
Пример: print("We have {} structures".format(e.size()))

get_structure(index: int) -> dict = {'energy': float, 'descr': str, 'xyz': np.array}
Возвращает данные о структуре с индексом index. Нумерация структур начинается с нуля, то есть 0 <= index <= e.size() - 1. Данные передаются в виде словаря с ключами 'energy' (энергия), 'descr' (описание) и 'xyz' (numpy матрица N на 3).
Пример: print("The first structure has energy {}".format(e.get_structure(0)['energy']))

get_atom_symbols() -> list of str
Возвращает лист из символов атомов.

Подсчет структур:

energy_count(e_cutoff: float, etype='kcal/mol' or 'hartree'/'a.u.') -> int
Возвращает количество структур в пределах данного e_cutoff от структуры с минимальной энергией (см. подробности в energy_filter).

distance_count(idxA: int, idxB: int, condition: function) -> None, condition(dist: float) -> bool
Возвращает количество структур, удовлетворяющих условию condition на длину связей между атомами idxA и idxB (см. подробности в distance_filter).

valence_count(idxA: int, idxB: int, idxC: int, condition: function) -> None, condition(angle: float) -> bool
Возвращает количество структур, удовлетворяющих условию condition на валентный угол между атомами idxA, idxB и idxC (см. подробности в valence_filter).

dihedral_count(idxA: int, idxB: int, idxC: int, idxD: int, condition: function) -> None, condition(angle: float) -> bool
Возвращает количество структур, удовлетворяющих условию condition на торсионный угол между атомами idxA, idxB, idxC и idxD (см. подробности в dihedral_filter).

TODO

1. RMSD filtering

2. Build connectivity graphs & search for VdW contacts

3. Complex filtering conditions (how to implement?)

4. Complex description templates (how to implement?)

TODO: include_from_dict

TODO: copy constructor, operator+