Задание 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 из предыдущего практикума.
$CONTRL COORD=CART 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
Только надо изменить информацию о типе входных координат: заменим COORD=CART на COORD=ZMT.
Проверим, работает ли наш файл-заготовка, т.е. запустим GAMESS. Выходной файл не содержит ошибок, а значит, переходим к следующему пункту.
Задание 3
Теперь создадим текстовый файл скрипта make_b.bash со следующим содержанием (скрипт сохраняем в формате UNIX во избежание всякого рода недоразумений):
#!/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. Запустим наш скрипт :
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
Непосредственно сам скрипт.
Результат работы скрипта - зависимость энергии молекулы от длины одной связи:
1.32986 -79.7339315289 1.34986 -79.7406173737 1.36986 -79.7462868578 1.38986 -79.7510331813 1.40986 -79.7549412921 1.42986 -79.7580886263 1.44986 -79.7605457819 1.46986 -79.7623771323 1.48986 -79.7636413832 1.50986 -79.7643920816 1.52986 -79.7646780790 1.54986 -79.7645439543 1.56986 -79.7640303992 1.58986 -79.7631745699 1.60986 -79.7620104053 1.62986 -79.7605689179 1.64986 -79.7588784564 1.66986 -79.7569649422 1.68986 -79.7548520851 1.70986 -79.7525615748 1.72986 -79.7501132582 |
Результат был сохранен в файле bond.
Задание 5
У нас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel. Но мы не ищем легких путей и сделаем это в gnuplot.Запустим Gnuplot:
gnuplot
Построим зависимость энергии от длины связи, просто введем : plot "bond".
У нас появился график с точками, похожими на параболу. Теперь нам надо найти коэффициенты в функции f(x)=a+k(x-b)^2, которые бы позволили наиболее близко описать наблюдаемую зависимость.
Сначала зададим функцию в развернутом виде, в строке 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)
Полученные значения:
Final set of parameters Asymptotic Standard Error ======================= ========================== a = -79.7652 +/- 0.0004522 (0.000567%) k = 0.563608 +/- 0.02335 (4.142%) b = 1.55432 +/- 0.002455 (0.1579%)Как видно, самая большая ошибка для k. Возможно, функция неточно совпадает с точками, поскольку в данном случае мы наблюдаем участок потенциала Мориса, а он имеет существенно более сложную зависимость.
Задание 6
Проделаем аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2.Все сделали точно так же, как было описано выше. Непосредственно скрипт.
Результат работы скрипта был сохранен в файле HCH.
В результате получили следующие коэффициенты:
Final set of parameters Asymptotic Standard Error ======================= ========================== 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.Итак, скрипт.
По полученным данным в gnuplot был построен график. В данном случае эта зависимость была аппроксимирована косинусом:
f(x)=k*cos(n*x)+a,
где k - константа жесткости, a - смещение, n - коэффициент растяжения/сжатия. Их начальные значения:
a = -80 n = 3 k = 0.05Результат работы скрипта был сохранен в файле d_3.
В результате получили следующие коэффициенты:
Final set of parameters Asymptotic Standard Error ======================= ========================== a = -79.7623 +/- 6.577e-07 (8.245e-07%) n = 3.00014 +/- 0.0002247 (0.007491%) k = 0.00234519 +/- 8.891e-07 (0.03791%)
А график выглядит вот как:
Мы имеем 3 минимума, так как крайние точки совпадают (углы -180 и 180).