Определение зависимости энергии молекулы этана от длины связи C-C
$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 с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не ее значением.
#!/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
Запустим скрипт:bash ./make_b.bash
В итоге получили 21 inp файл с разными значениями для переменной сс.
gms b_${i}.inp 1 > b_${i}.log
Для извлечения значения энергий из log файла закомментируем в скрипте запуск Gamess. Затем добавим вызов awk:
echo -n "$nb " (pаспечатаем переменную nb и несколько пробелов без переноса строки);
awk '/TOTAL ENERGY =/{print $4}' b_${i}.log (ищем строчку с TOTAL ENERGY и печатаем 4ое поле считая что поля разделены пробелами).
После запуска скрипта перенаправим в файл две колонки цифр (зависимость энергии от длины связи):
bash ./make_b.bash > bond
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
Сохраним значения коэффициентов:
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%)
correlation matrix of the fit parameters:
a k b
a 1.000
k -0.724 1.000
b 0.171 -0.413 1.000
Построим графики функции и значений энергии из Gamess:
plot "bond", f(x)
Cамая большая ошибка для k. Возможно, из-за этого полученная функция и точки совпадают не абсолютно.
Определение зависимости энергии молекулы этана от значения валентного угла HCH
Коэффициенты f(x), подогнанные gnuplot (сама функция и стартовые значения коэффициентов аналогично предыдущему заданию):
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%)
Определение зависимости энергии молекулы этана от значения торсионного угла d3
По полученным данным в gnuplot был построен график. Он похож на синусоиду, поэтому эта зависимость была аппроксимирована следующей функцией:
f(x)=k*cos(n*x)-a
где k - константа жесткости, a - смещение, n - коэффициент растяжения/сжатия. Их начальные значения:
a = -80
n = 3
k = 1
После подгонки коэффициенты приобрели следующий вид:
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%)
График, полученный в gnuplot:
Мы наблюдаем 3 минимума, так как крайние точки совпадают (угол -180 равне углу 180).
Определение зависимости энергии молекулы этана от длины связи C-C при увеличении шага до 0.1 ангстрема
полученые значения
Значения коэффициентов:
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = -79.7647 +/- 5.444e-05 (6.825e-05%)
k = 0.53427 +/- 0.01117 (2.091%)
b = 1.53859 +/- 0.0005931 (0.03855%)
correlation matrix of the fit parameters:
a k b
a 1.000
k -0.734 1.000
b 0.126 -0.308 1.000 
Аппроксимация получилась ровно такая же (и функция, и значения параметров), как и в самом первом задании. Но связь получилась короче и энергия повыше (в первом задании она была в диапазоне [-79,77...-79,73], здесь же [-79,765...-79,758]).