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

Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчетов.
1. Предоставлена оптимизированная структура этана в виде 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. Bash это интерпретатор командной строки который автоматически запускается при присоединении к kodomo.



2. Составим файл-заготовку для размножения et.inp. Для этого к координатам добавим шапку для dft из предыдущего практикума (только надо изменить информацию о типе входных координат: заменим COORD=CART на COORD=ZMT):
   

   $CONTRL COORD=ZMT UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY $END
   $BASIS  GBASIS=N31 NGAUSS=6
   POLAR=POPN31 NDFUNC=1 $END
   $GUESS  GUESS=HUCKEL $END
   $system mwords=2 $end
   $DATA
           




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

   #!/bin/bash
   ### делаем цикл от -10 до 10 ##### 
   
   for i in {-10..10}; do 
   
   #### нам надо рассчитать новую длину связи #####
   #### с шагом 0.02 ангстрема,               #####
   #### воспользуемся калькулятором bc        #####
   #### и результат поместим в переменную nb  #####
   
        nb=$(echo "scale=1.52986; 1.001 + $i/50" | bc -l)
   
   #### пролистаем файл et.inp и заменим указание переменной ###
   #### на новое значение и пере направим результат в файл   ###
   
        sed "s/cc=1.001/cc=$nb/" et.inp  > b_${i}.inp
   
   done
    

   где 1.52986 - стартовая длина изменяемой связи. Запустим скрипт:
   
   bash ./make_b.bash
   
   В результате получаем 21 inp файл, в каждом - разное значение для переменной сс.


4. Теперь запустим расчет для этих файлов. Для этого перед строчкой с

   done
   

   вставим запуск Gamess:

   gms b_${i}.inp 1 > b_${i}.log
   

   Запустим скрипт. Теперь надо извлечь значение энергии из log файла. Удобно воспользоваться awk. Сначала в скрипте комментируем запуск Gamess поставив в начало строчки c gms # . Добавим после за комментированной строчки вызов awk, при этом мы ищем строчку с TOTAL ENERGY и печатаем 4ое поле считая что поля разделены пробелами:

   awk '/TOTAL ENERGY =/{print $4}'  b_${i}.log 
   
   

   Запускаем скрипт. Получаем 21 значение энергии.
   Теперь удобно было бы выводить и значение длины связи, для этого добавим перед вызовом awk распечатку переменой nb. Распечатаем переменную и несколько пробелов без переноса строки:

   
   echo -n "$nb"
   
       

   Теперь запустиМ скрипт и перенаправляем поток в файл : bash ./make_b.bash > bond.


5. Таким образом, мы имеем зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel. Предлагается сделать это в gnuplot.
   Запустим Gnuplot:

   
   gnuplot
              

   Построим зависимость энергии от длины связи Появляеться график с точками похожими на параболу. Теперь надо найти коэффициенты в функции 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 
   

   Значения коэффициентов:

   a               = -79.7652         +/- 0.0004522    (0.000567%)
   k               = 0.563608         +/- 0.02335      (4.142%)
   b               = 1.55432          +/- 0.002455     (0.1579%)
   

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


6. Проделаем аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2. Построим графики функции и значений энергии из Gamess:
   
   
   
   Значения коэффициентов:

   a               = -79.7647         +/- 1.21e-08     (1.517e-08%)
   k               = 3.56076e-05      +/- 6.229e-09    (0.01749%)
   b               = 111.38           +/- 9.954e-05    (8.937e-05%)
   



7. Проделаем аналогичные операции для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12.
   
   
   
   Значения коэффициентов:

   	
   a               = 0.00234519       +/- 8.891e-07    (0.03791%)
   k               = 3.00014          +/- 0.0002247    (0.007491%)
   b               = -79.7623         +/- 6.577e-07    (8.245e-07%)
   

   Видим, что число минимумов не меньше 3-х.