$BASIS GBASIS=N31 NGAUSS=6
POLAR=POPN31 NDFUNC=1 $END
$GUESS GUESS=HUCKEL $END
$system mwords=2 $end
$DATA
Чтобы проверить, работает ли файл, я запустила Gamess.
3.
Создаем файл скрипта make_b.bash, как описано в задании, поменяв значение стартовой длины изменяемой связи так, чтобы она соответствовала файлу et.inp:
#!/bin/bash
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)
sed "s/cc=1.001/cc=$nb/" et.inp > b_${i}.inp
done
Запускаем скрипт: bash ./make_b.bash
На выходе получаем 21 inp файл, при этом в каждом разное значение для переменной cc.
Затем необходимо запустить рассчет для полученных файлов. Для этого приписываем к файлу скрипта перед строкой done следующее:
gms b_${i}.inp 1 > b_${i}.log
Чтобы извлечь значения инергии из log файлов, используем команду awk. Перед строкой gms ставим знак # (комментируем), а после вводим:
awk '/TOTAL ENERGY =/{print $4}' b_${i}.log
Чтобы вывести значения длин связей, ставим следующее перед командой awk:
echo -n "$nb "
Запускаем скрипт и пишем результат в файл:
bash ./make_b_awk.bash > bond
Результаты находятся в файле bond
5.
Далее необходимо построить зависимость энергии молекулы от длины одной связи. Для начала запустила gnuplot и ввела команду: plot "bond"
Рис. 1. Длина связи C-C.
На рисунке 1 видим график, похожий на параболу. Ищем коэффициенты в функции f(x)=a+k(x-b)^2, которые бы позволили наиболее близко описать наблюдаемую зависимость. Задаем функцию в развернутом виде: 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.
Полученные коэффициенты лежат в файле fit1.log
Рис. 2. Длина связи C-C.
Как можно видеть на рисунке 2, функция не совсем совпадает. Можно предположить, что в данном случае нужно использовать другую, более сложную функцию зависимости энергии связи от ее длины.
6.
Необходимо проделать те же операции для валентного угла HCH, при этом его значения должны изменяться от 109.2 до 113.2.
Результаты находятся в файле bond2
Рис. 3. Валентный угол HCH.
7.
Повторяем процедуру для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12. Не проводим подгонку.
Результаты находятся в файле bond3
Рис. 4. Торсионный угол d3.
Важно понимать, что -180 и +180 - это один и тот же угол, а значит, мы видим три минимума и три максимума.
©Melnichuk Anastasia