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

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

Необходимые сведения о работе с GAMESS см. здесь.
Сведения о работе с Gnuplot см. здесь.
Введения о скриптовании  в Bash здесь.
Ведения о awk здесь.
Вся работа по расчётам будет проходить на kodomo через терминал putty, а для работы с графическим выводом Gnuplot понадобится Xming.
Пример подобной работы представлен в статье о разработке поля gaff.

  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. Cоставим файл-заготовку для размножения. 
    Для этого к координатам добавим шапку:

 $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 
Проверим работает ли файл-заготовка (запустим GAMESS, и проверим выходной файл на наличие ошибок).

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

    #!/bin/bash
### делаем цикл от -10 до 10 ##### 
 
for i in {-10..10}; do 
 
#### нам надо рассчитать новую длину связи #####
#### с шагом 0.02 ангстрема,               #####
#### воспользуемся калькулятором bc        #####
#### и результат поместим в переменную nb  #####
 
     nb=$(echo "scale=5; 1.001 + $i/50" | bc -l)
 
#### пролистаем файл et.inp и заменим указание переменной ###
#### на новое значение и пере направим результат в файл   ###
 
     sed "s/cc=1.001/cc=$nb/" et.inp  > b_${i}.inp
 
done
    В скрипте стартовая длина изменяемой связи соответствовует файлу et.inp. Запустим скрипт : 
     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. У нас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в gnuplot.
    Для этого запустим Xming->XLaunch. Выберем тип расположения окон, удобно использовать Multiple windows. Next. Выберем Start Program. Run Remote-> Putty. Дальше всё как обычно: kodomo, username.
    Перейдем в рабочую директорию. Запустим Gnuplot:

    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
    Сохраните значения коэффициентов. Постройте графики функции и значений энергии из Gamess. 
    plot "bond", f(x)

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

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

    Параметры подгонки. Полученная апроксимация хорошая - все точки лежат прямо на кривой.

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

    Подгонку не проводили. Число минимумов функции как минимум 3.

  8. Увеличим шаг до 0.1 ангстрема при расчёте связи сс. Для этого запустим скрипт и перенаправим поток в файл. Построим зависимость значений энергии из Gamess:

    Параметры подгонки. Как видно наблюдаемую зависимость можно аппроксимировать экспонентой.

© Пискунова Юлия 2011