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

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

Имеем оптимизированную структуру этана в виде z-matrix:
```
 $DATA
eth
C1
 C   
 C      1   cc 
 H      2   ch   1   cchv 
 H      2   ch   1   cch   3   d1 0
 H      2   ch   1   cch   3   d2 0
 H      1   ch   2   cch   3   d3 0
 H      1   ch   2   cch   5   d3 0
 H      1   chv  2   cch   4   d3 0

cc=1.52986
ch=1.08439
chv=1.08439
cch=111.200
cchv=111.200
d1=120
d2=-120
d3=180
 $END
```
Вместо значений длин и углов связей стоят переменные. Необходимо создать порядка 20 разных файлов для расчёта энергии в Gamess с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не её значением. Автоматизируем процесс с помощью скрипта на bash.
Создадим файл-заготовку для размножения et.inp. Для этого к координатам добавим шапку для DFT из предыдущего практикума. Нужно изменить информацию о типе входных координат: заменим COORD=CART на COORD=ZMT.
Проверим, работает ли файл-заготовка. Для этого запустим GAMESS. В результате получаем файл et.log без ошибок. Значит, можно двигаться дальше.

Теперь создадим текстовый файл скрипта make_b.bash со следующим содержанием:

#!/bin/bash
### делаем цикл от -10 до 10 ##### 

for i in {-10..10}; do 

#### нам надо рассчитать новую длину связи #####
#### с шагом 0.02 ангстрема,               #####
#### воспользуемся калькулятором bc        #####
#### и результат поместим в переменную nb  #####

     nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)

#### пролистаем файл et.inp и заменим указание переменной ###
#### на новое значение и пере направим результат в файл   ###

     sed "s/cc=1.52986/cc=$nb/" et.inp  > b_${i}.inp

done

В нем 1.52986 - стартовая длина изменяемой связи (взята из et.inp).
Запустим скрипт :

bash ./make_b.bash

В результате получаем 21 inp файл (в каждом - разное значение для переменной сс).

Запустим расчёт для этих файлов.
Для этого перед строчкой
```
done
```
вставим запуск Gamess:
```
gms b_${i}.inp 1 > b_${i}.log
```
Запустим скрипт.

Теперь необходимо извлечь значение энергии из log файла. Для этого удобно воспользоваться awk.
Сначала в скрипте прокомментируем запуск Gamess, поставив в начало строчки
```
     gms b_${i}.inp 1 > b_${i}.log
```
####.

Добавим после этой строчки вызов awk (мы ищем строчку с TOTAL ENERGY и печатаем четвертое поле, считая, что поля разделены пробелами):
```
awk '/TOTAL ENERGY =/{print $4}'  b_${i}.log
```
Теперь удобно было бы выводить и значение длины связи. Для этого добавим перед вызовом awk распечатку переменой nb. Распечатаем переменную и несколько пробелов без переноса строки:
```
echo -n "$nb    "
```
Наконец, перенаправим поток вывода скрипта в файл bond:
```
bash ./make_b.bash > bond
```
В результате имеем файл bond с двумя столбцами (с длинами связи и соответствующими значениями энергии).
Получаем зависимость энергии молекулы от длины одной связи. Построим ее в gnuplot. Для этого просто введем:
```
plot "bond"
```
Получаем график с точками, похожими на параболу. Теперь надо найти коэффициенты в функции 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
```
Проведём подгонку коэффициентов под имеющиеся точки в файле bond:
```
fit f(x) "bond" via a,k,b
```
Значения коэффициентов сохраним в файле koeff.
Построим графики функции и значений энергии из Gamess:
```
plot "bond", f(x)
```
Изображение графика представлено ниже:
Проделаем аналогичные операции для валентного угла HCH (cch). Будем изменять его значения от 109.2 до 113.2 с шагом 0.2. В результате получим файл angle со значениями угла и соответствующими им энергиями. Зависимость энергии от угла была аппроксимирована параболой со значениями коэффициентов koeff2. График зависимости представлен ниже:
Проделаем аналогичные операции для торсионного угла d3. Будем изменять его значения от -180 до 180 с шагом 12. В результате получим файл tors со значениями угла и соответствующими им энергиями. Зависимость энергии от угла представлена ниже:

Как видно, эта функция имеет 3 минимума (если считать, что угол -180 равен углу 180).
Рассчитаем связь, увеличив шаг до 0.1 ангстрема. В результате получим файл bond2 со значениями связи и соответствующими им энергиями. Зависимость энергии от угла была аппроксимирована немного модифицированным потенциалом Морзе (f(x)=a*(exp((-2)*b*(x-c))-2*exp((-1)*b*(x-c))+exp(d*x))) со значениями коэффициентов koeff4. График зависимости представлен ниже: