Вычисление параметров для молекулярной механики

Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчетов.

У нас имеется оптимизированная структура этана, представленная в виде 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).