Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.
У нас имеется оптимизированная структура этана, представленная в виде z-матрицы:
Вместо значений длин и углов связей стоят переменные.
Наша цель состоит в том, что бы создать порядка 20 разных файлов для расчёта энергии в Gamess с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не её значением.
Попробуем автоматизировать процесс с помощью скрипта на bash.
1. Создадим файл-заготовку для размножения. et.inp
Для этого к координатам добавьте шапку для dft из предыдущего практикума:
Не забудем изменить информацию о типе входных координат: заменим COORD=CART на COORD=ZMT!!!
Запустим GAMESS:
gms et.inp 1 >& et.log
Выходной файл: et.log ошибок не содержит. Поэтому приступаем к следующему шагу.
2. Составим скрипт на bash, чтобы сгенерировать 21-inp файл с разными длинами связи с-с.
Теперь добавим строку, для того, чтобы провести расчет энергии с помощью 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
Полученный скрипт: make_b.bash .
Для запуска скрипта используется команда:
bash ./make_b.bash
Подобные операции проводим для валентного угла HCH (make_b.bash), его значения должны изменяться от 109.2 до 113.2 , а также для для торсионного угла d3 (make_b.bash), его значения должны изменяться от -180 до 180 c шагом 12.
3. У вас есть зависимость энергии молекулы от длины одной связи С-С : bond.
(также зависимость энергии от изменений валентного угла HCH (bondНСН) и торсионного угла d3 (d3bond))
Построим эту зависимость с помощью gnuplot:
Запускаем Xming для Putty и присоединяемся к kodomo. Переходим в рабочую директорию и запускаем Gnuplot:
gnuplot
Построим зависимость энергии от длины связи:
plot "bond"
Появился график с точками, похожими на параболу. Теперь нам надо найти коэффициенты в функции 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
Сохраните значения коэффициентов. Построим графики функции и значений энергии из Gamess:
plot "bond", f(x)
файл с коэффициентами: kocc
Функция неточно совпадает с точками, так как в данном случае мы наблюдаем участок потенциала Мориса, имеющий более сложную зависимость.
Аналогично, график был построен для манипуляций с валентными углом НСН:
файл с коэффициентами: koНСН
Данная зависимость хорошо апроксимируется параболической функцией.
Подобная операция, проведённая для торсионного угла d3, дала следующий результат:
файл с коэффициентами: kod3
(Апроксимация проводилась функцией: f(x)=a*cos(k*x/180*pi) + b )
Мы имеем 3 минимума, так как крайние точки совпадают (углы -180 и 180).