запускаем 'ipython notebook' и 'pymol -R' из рабочей директории.
Форматирование задаётся Markdown: http://daringfireball.net/projects/markdown/basics

сведения о Морас можно получить здесь: http://openmopac.net/manual/index.html

Работа проводится с помощью putty и обработки pymol. Команды для bash выделены табуляцией.

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


Задание 1. Оптимизация 3D структур нафталина и азулена.


сперва установим переменные: 

export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
In [1]: 
#Для работы с Pymol: прокси и загрузка картинок:
from xmlrpclib import ServerProxy
from IPython.display import Image
import time

cmd = ServerProxy(uri="http://localhost:9123/RPC2")

фишка в том, что я запускаю ipython c Windows, в виду того, что порты на кодомо закрыты, это не предосудительно. У меня нет Bash, потому я просто копирую строки из скрипта ниже в putty, и они выполняются по мере готовности. Тем не менее, тэг "%%bash" я поставлю, вдруг из класса будет запускаться.

Нам дана z-матрица этана, на основании которой надо сделать 21 файл с варьируемыми параметрами (начнем с длин связи).

сделаем файл из заголовка с прошлого практикума, где заменим COORD=CART на COORD=ZMT, после чего пропишем полностью нормальную z-матрицу этана.

In [2]: 
%%bash
BASE="\$CONTRL COORD=ZMT UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY \$END\n\
 \$BASIS  GBASIS=N31 NGAUSS=6\n\
  POLAR=POPN31 NDFUNC=1 \$END\n\
 \$GUESS  GUESS=HUCKEL \$END\n\
 \$system mwords=2 \$end\n\
 \$DATA\n\
 eth\n\
 C1\n\
  C   \n\
  C      1   cc \n\
  H      2   ch   1   cchv \n\
  H      2   ch   1   cch   3   d1 0\n\
  H      2   ch   1   cch   3   d2 0\n\
  H      1   ch   2   cch   3   d3 0\n\
  H      1   ch   2   cch   5   d3 0\n\
  H      1   chv  2   cch   4   d3 0\n\
  \n\
 cc=1.52986\n\
 ch=1.08439\n\
 chv=1.08439\n\
 cch=111.200\n\
 cchv=111.200\n\
 d1=120\n\
 d2=-120\n\
 d3=180\n\
  \$END"

#запишем это в файл

echo -en " " > et.inp
echo -e $BASE >> et.inp

#и проверим, работает ли наша химера

gms et.inp  1 >& et.log

#у меня работает, файл по сравнению с прошлыми log-файлами, запущенными с другими параметрами получился короче, считается быстрее, 
#наверное, дело в меньшем числе атомов, все-таки это этан.

#ниже вы видите, что bash у меня нет, см. комментарий вверху страницы курсивом.

Couldn't find program: u'bash'

далее сделаем скрипт на bash , чтобы создать 21 файл с варьируемой длиной связи СС по образцу на https://kodomo.fbb.msu.ru/wiki/2011/8/task5 , изменив длины связей на подходящие для этана в соответствии с реальной z-матрицей (выше).

В коде скрипта строки до закомменченного запуска GAMESS отвечают за создание 21 файла с разной длиной связи сс, активированная строка с GAMESS в цикле запускает создание файлов .log с вычислениями, а строки под строкой GAMESS используют эти файлы для того, чтобы вывесли длину связи СС и через несколько пробелов общую энергию. Эту выдачу мы записываем в файл командой

bash ./make_b.bash > bond.txt


Проблема с windows возникает, если писать скрипт и сохранять файл не через системы linux, тогда переносы строк содержат еще символ , который не видно, и надо скрипт после сохранения прогнать например через 

sed -e 's/\r$//' inputfile > outputfile

альтернативно люди предлагают использовать dos2unix.

Теперь используем Xming->XLaunch -> Multiple windows -> Start Program -> Putty -> kodomo, username -> Переходим в рабочую директорию, Запускаем Gnuplot:

gnuplot

Строим зависимость энергии от длины связи:

plot "bond.txt"

Появились точки, похожие на параболу. Теперь нам надо найти коэффициенты в функции f(x)=a+k(x-b)^2 которые бы позволили наиболее близко описать наблюдаемую зависимость. Для этого воспользуемся возможностями Gnuplot. Сначала зададим функцию в развернутом виде, в строке gnuplot введём:

f(x)=a + k*x*x - 2*k*x*b + k*b*b

И зададим стартовые значения коэффициентов:

a=-80
k=1
b=1.5

Настроим gnuplot для сохранения изображения в файл (название файла и параметры задаются, файл сохраняется в рабочую директорию).

set terminal png size 1000,700 enhanced font "Helvetica,16"
set output 'output.png'

Проведём подгонку коэффициентов под имеющиеся точки в файле bond:

fit f(x) "bond.txt" via a,k,b
В результате получаем параметры 
a = -79.7652
k = 0.5636
b = 1.5543

Для совместного изображения линии и графика используем

plot "bond.txt", f(x)

ссылка на картинку ниже

In [10]: 
from IPython.display import display_html
Image(url='http://kodomo.cmm.msu.su/~aleshin_vasily/projects/Term8/images/output_cc.png')
Out[10]:

Рис. 1. Зависимость общей энергии этана от длины связи СС (значения отмечены точками). Зеленым показана аппроксимирующая квадратичная функция (коэффициенты приведены выше).

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

Чтобы проделать то же с углом НСН заменим в скрипте переменную сс на chv (см. исходную матрицу выше). Исходно все углы HCH принимают значение 111.2, мы же проварьируем их от 109.2 до 113.2. Возьмем шаг 0.2 градуса.

Переделаем и используем новый скрипт для создания 21 файла с разными углами и их обсчета в GAMESS.

bash ./make_ch.bash

Далее лишь закомментируем строку с использование GAMESS в этом скрипте и откроем две строки под ней для того, чтобы выписать значения энергии из файлов.

bash ./make_ch.bash > angle.txt

Новый файл со значениями одного угла НСН и энергией этана визуализируем:

Снова воспользуемся Xming->XLaunch для запуска Gnuplot.

gnuplot

График снова похож на параболу, потому снова можно аппроксимировать той же функцией, но начальные параметры слегка изменим:

f(x)=a + k*x*x - 2*k*x*b + k*b*b  
a=-80
k=1

b=111.2

fit f(x) "bond.txt" via a,k,b
В результате получаем параметры 
a = -79.7647
k = 3.56076e-05
b = 111.38

set terminal png size 1000,700 enhanced font "Helvetica,16"
set output 'output2.png'

plot "angle.txt", f(x)
In [5]: 
from IPython.display import display_html
Image(url='http://kodomo.cmm.msu.su/~aleshin_vasily/projects/Term8/images/output2.png')
Out[5]:

Можно заметить странную вещь - минимум энергии при расчете пришелся на значение угла 111.38, а не 111.2, как это приведено в изначальной молекуле. Возможно, дело в том, что мы использовали лишь один из углов, но скорее дело в разных использованных полях при расчете энергии GAMESS.

Напоследок варьируем один из торсионных углов - d3 от -180 до 180 с шагом 12 используя скрипт, где сперва используем строку с запуском GAMESS:

bash ./make-tor.bash

затем закомменченные строки вместо GAMESS:

bash ./make-tor.bash > angle_tor.txt

И получаем файл с углами и энергиями.

Используем Х-сервер и gnuplot:

set terminal png size 1000,700 enhanced font "Helvetica,16"
set output 'output-tor.png'

plot "angle-tor.txt"
In [10]: 
from IPython.display import display_html
Image(url='http://kodomo.cmm.msu.su/~aleshin_vasily/projects/Term8/images/output_tr.png')
Out[10]:
Как и стоило ожидать, энергия молекулы зависит от торсионных углов как периодическая функция с периодом 120 градусов, что хорошо заметно на рисунке выше.

в финале конвертируем файл ipython в html:

ipython nbconvert --to html pract4.ipynb