Занятие 4.
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
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.
2.Составим файл-заготовку для размножения, пусть имя файла будет et.inp. Для этого к координатам добавьте шапку для dft из предыдущего практикума. Только надо изменить информацию о типе входных координат: замените COORD=CART на COORD=ZMT.
3.Теперь создадим текстовый файл скрипта 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
Поправим скрипт так, что бы стартовая длина изменяемой связи соответствовала файлу et.inp. make_b.bash Запустим скрипт : Это и будут входные файлы для оптимизации геометрии средствами GAMESS.bash ./make_b.bashПоявился 21 inp файл и в каждом разное значение для переменной сс.
2.Запустим расчёт для этих файлов, для этого перед строчкой с 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 "Теперь запустим скрипт и если результат подходит, т.е. есть две колонки цифр (bond), то можно пере направить поток в файл :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Сохраним значения коэффициентов.коэффициенты Большая ошибка для k, из-за этого функция и точки совпадают не точно. Построим графики функции и значений энергии из Gamess.plot "bond", f(x)
6.Проделаем аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2. make_b_2.bash углы_энергии Сохраним полученные коэффициенты в отчёт.коэффициенты Зависимость энергии от угла аппроксимировали параболой:
7.Проделаем аналогичные операции для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12. make_b_3.bash углы_энергии Сохраним полученные коэффициенты в отчёт. коэффициенты Зависимость энергии от угла аппроксимировали косинусом f(x)=k*cos(n*x)+a (k=0,05, n=3, a=-80) Функция имеет три минимума:
