Учебный сайт Фоменко Елены
| Главная | Семестры | Проекты | Заметки |
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
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
2. Создадим порядка 20 разных файлов для расчёта энергии в Gamess с разными значениями по длине одной из связей.
Попробуем автоматизировать процесс с помощью скрипта на bash.
Составляем файл-заготовку для размножения: et.inp.
Для этого к координатам добавили шапку для dft из предыдущего практикума и изменили информацию
о типе входных координат, COORD=CART на COORD=ZMT:
$CONTRL COORD=ZMT 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
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
Запускаем GAMESS. Выходной файл не содержит ошибок, так что переходим к следующему пункту.
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
Запускаем скрипт:
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
Получаем файл bond с двумя колонками чисел.
5. Таким образом, у нас есть зависимость энергии молекулы от длины одной связи. Строим эту зависимость в 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
Проведём подгонку коэффициентов под имеющиеся точки:
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%)
Построим графики функции и значений энергии из Gamess (рис.1):
plot "bond", f(x)
Рис.1 Зависимость энергии молекулы этана от длины связи C-C и результат аппроксимации функцией f(x).
6. Проделаем аналогичные операции для валентного угла HCH (переменная cchv), его значения должны изменяться от 109.2 до 113.2 (скрипт: make_b.bash2).
Рис.2 Зависимость энергии молекулы от величины валентного угла HCH и результат ее аппроксимации функцией f(x).
С прежними начальными коэффициентами (рис.2) аппроксимация не прошла, при их изменении (a=-70, k= 10, b=110) все получилось, причем точнее, чем в предыдущий раз:
Рис.3 Зависимость энергии от величины валентного угла HCH и результат аппроксимации зависимости
при измененных начальных коэффциентах.
Полученные итоговые коэффициенты:
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = -79.7647 +/- 1.21e-008 (1.517e-008%)
k = 3.56076e-005 +/- 6.229e-009 (0.01749%)
b = 111.38
7. Сделаем то же самое для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12 (скрипт: make_b.bash3).
Рис.4 Зависимость энергии от величины торсионного угла d3.
Количество минимумов функции - 3 (180 и -180 аналогичны).