Учебная страница курса биоинформатики,
год поступления 2019

Семестры Студенты Преподаватели

• 2024
• 2023
• 2022
• 2021
• 2020
• 2018
• 2017
• 2016

• ...

Вычисление параметров для молекулярной механики

• Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов. Пример подобной работы представлен в статье о разработке поля gaff.

   1 Если задание будет сделано с помощью ASE то двойной бал. https://wiki.fysik.dtu.dk/ase/dev/index.html

Для реализации задания на Google Colab:

   1 !pip install -q condacolab
   2 import condacolab
   3 condacolab.install()
   4 condacolab.check()

   1 ! mamba install -c anaconda intel-openmp
   2 ! mamba install -c psi4 psi4

   1 import psi4
   2 import numpy as np
   3 psi4.core.set_output_file('output.dat')

Расчёт орбиталей в Psi4: информация об api https://psicode.org/psi4manual/1.4.0/api/psi4.core.cubeproperties

Основное задание

1. Вам предоставлена оптимизированная структура этана в виде z-matrix :

   1 inp = '''!HF RHF 6-31G
   2 * int 0 1
   3 C 0 0 0 0 0 0 
   4 C 1 0 0 1.52986 0 0 
   5 H 1 2 0 1.08439 111.200 0
   6 H 1 2 3 1.08439 111.200 120
   7 H 1 2 3 1.08439 111.200 -120
   8 H 2 1 3 1.08439 111.200 180
   9 H 2 1 5 1.08439 111.200 120
  10 H 2 1 5 1.08439 111.200 -120
  11 *
  12 '''

Как вы видите, вместо значений длин и углов связей можно вставить значение. Наша цель состоит в том, что бы рассчитать 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей (C-C).

2. Конечно можно сделать эти файлы вручную, но ожидается, что вы сделаете это средствами python .

3. Предлагаю следующую функцию по запуску Psi4

   1 def run_psi4(inp):
   2     m = psi4.geometry(g)
   3     psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
   4     psi4.energy('scf/cc-pvtz', molecule = m )
   5 
   6     return ....

• Можно посмотреть

import py3Dmol
view = py3Dmol.view()
view.addModel(m.save_string_xyz_file(),'xyz')

4. У вас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel. Вам предлагается сделать это в matplotlib в Jupyter Notebook.
   

Для matplotlib:

   1 %matplotlib inline
   2 import numpy as np
   3 import matplotlib.pyplot as plt
   4 from scipy import optimize 

   1 # fake x array , replace with real one
   2 x_o=np.arange(1,2,0.1)
   3 # fake y array, replace with energies
   4 y_o=-70+2*(x_o -1.58)**2 
   5 
   6 #function is  f(x)=k(b-x)^2 + a
   7 fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
   8 errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
   9 
  10 p0 = [1,1, -79] # Initial guess for the parameters
  11 p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
  12 print "Optimized params:", p1
  13 
  14 #Plot it
  15 plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
  16 plt.xlim(1,2)
  17 plt.show()

5. Проделайте аналогичные операции для валентного угла HCС, его значения должны изменяться от 109.2 до 113.2. Сохраните полученные коэффициенты в отчёт. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчёта. Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.
6. Проделайте аналогичные операции для торсионного угла CC, его значения должны изменяться от -180 до 180 c шагом 12. Сохраните изображение и укажите в отчёте количество минимумов функции.
7. * Увеличьте шаг до 0.1 ангстрема при расчёте связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.
8. * Установите psiresp и рассчитайте точечные заряды для этана.

mamba install -c conda-forge psiresp