|
Занятие 4. Вычисление параметров для молекулярной механики.
Определение зависимости энергии молекулы этана от длины связи C-C.
Оптимизированная структура этана в виде 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 с разными значениями по длине одной из связей.
Файл-заготовка для размножения: et.inp.
Cкрипт make_b.bash получает 21 inp файла и в каждом разное значение для переменной сс:
#!/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
Затем с помощью этого же скрипта для полученных inp файлов был произведен расчет программой Gamess. Для этого в скрипт была добавлена строчка:
gms b_${i}.inp 1 > b_${i}.log
Запустите скрипт и через какое-то время расчет закончится. Теперь нам надо извлечь значение энергии из
log файла. Удобно воспользоваться awk.
Значения энергии из log файлов были извлечены с помощью awk. Для этого снова был использован скрипт 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
# gms b_${i}.inp 1 > b_${i}.log
echo -n "$nb "
awk '/TOTAL ENERGY =/{print $4}' b_${i}.log
done
|
Результат работы скрипта - зависимость энергии молекулы от длины одной связи:
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.
С помощью программы gnuplot была изображена зависимость энергии от длины связи по данным из файла bond.
Т.к. полученный график похож на параболу, данная зависимость была аппроксимирована функцией f(x)=a+k(x-b)^2.
В Gnuplot была совершена подгонка коэффициентов под имеющиеся точки в файле bond (сначала была задана функция, затем
стартовые значения коэффициентов и только потом сама подгонка):
f(x)=a + k*x*x - 2*k*x*b + k*b*b
a=-80
k=1
b=1.5
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%)
Как видно, самая большая ошибка для k. Возможно, из-за этого полученная функция и точки совпадают не абсолютно.
Определение зависимости энергии молекулы этана от значения валентного угла HCH.
Для валентного угла HCH были проведены аналогичные операции, причем его занчения изменялись от 109.2 до 113.2
(смотреть скрипт):
109.20000 -79.7645100655
109.40000 -79.7645396319
109.60000 -79.7645663700
109.80000 -79.7645902763
110.00000 -79.7646113475
110.20000 -79.7646295805
110.40000 -79.7646449722
110.60000 -79.7646575202
110.80000 -79.7646672219
111.00000 -79.7646740754
111.20000 -79.7646780790
111.40000 -79.7646792311
111.60000 -79.7646775309
111.80000 -79.7646729776
112.00000 -79.7646655711
112.20000 -79.7646553114
112.40000 -79.7646421992
112.60000 -79.7646262354
112.80000 -79.7646074214
113.00000 -79.7645857589
113.20000 -79.7645612501
|
График, полученный в gnuplot:
Коэффициенты f(x), подогнанные gnuplot (сама функция и стартовые значения коэффициентов аналогично пункту 1):
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%)
Сохраните полученные коэффициенты в отчет. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчета.
Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.
Определение зависимости энергии молекулы этана от значения торсионного угла d3.
Для торсионного угла в3 были проведены аналогичные п.1. и п.2 операции, причем его занчения изменялись от 180 до -180
(смотреть скрипт):
180 -79.7646780790
168 -79.7642336928
156 -79.7630628377
144 -79.7616148156
132 -79.7604407798
120 -79.7599827599
108 -79.7604407798
96 -79.7616148156
84 -79.7630628377
72 -79.7642336928
60 -79.7646780790
48 -79.7642336928
36 -79.7630628377
24 -79.7616148156
12 -79.7604407798
0 -79.7599827599
-12 -79.7604407798
-24 -79.7616148156
-36 -79.7630628377
-48 -79.7642336928
-60 -79.7646780790
-72 -79.7642336928
-84 -79.7630628377
-96 -79.7616148156
-108 -79.7604407798
-120 -79.7599827599
-132 -79.7604407798
-144 -79.7616148156
-156 -79.7630628377
-168 -79.7642336928
-180 -79.7646780790
|
По полученным данным в 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:
Как видно, у функции f(x) на отрезке от -180 до 180 три минимума.
|
|
|