Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.cmm.msu.ru/~anuta_al/MolDin4.html
Дата изменения: Mon May 30 16:17:17 2011
Дата индексирования: Tue Oct 2 01:31:44 2012
Кодировка: Windows-1251
Вычисление параметров для молекулярной механики

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



Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчетов.
  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
    
    Вместо значений длин и углов связей стоят переменные. Цель состоит в том, чтобы создать порядка 20 разных файлов для расчета энергии в Gamess с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не ее значением.
    Автоматизируем процесс с помощью скрипта на bash.
  2. Составим файл-заготовку для размножения, пусть имя файла будет et.inp.
    Для этого к координатам добавим шапку для dft из предыдущего практикума. Только надо изменить информацию о типе входных координат: заменим COORD=CART на COORD=ZMT.
    Проверим, работает ли файл-заготовка: выходной файл et.log не содержит ошибок, переходим к следующему пункту.

  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.001/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, поставив в начало строчки
         gms b_${i}.inp 1 > b_${i}.log
    ####. 
    
    Добавим после этой строчки вызов awk (мы ищем строчку с TOTAL ENERGY и печатаем четвертое поле, считая, что поля разделены пробелами):
    awk '/TOTAL ENERGY =/{print $4}'  b_${i}.log
    Теперь удобно было бы выводить и значение длины связи. Для этого добавим перед вызовом awk распечатку переменой nb. Распечатаем переменную и несколько пробелов без переноса строки:
    echo -n "$nb    "
    Перенаправим поток вывода скрипта в файл bond:
    bash ./make_b.bash > bond
    В результате имеем файл bond с двумя столбцами (с длинами связи и соответствующими значениями энергии).

  5. Получаем зависимость энергии молекулы от длины одной связи. Построим ее в 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
    Значения коэффициентов сохраним в файле koeff.txt.


    Построим графики функции и значений энергии из Gamess:
    plot "bond", f(x)
    Изображение графика:





  6. Проделаем аналогичные операции для валентного угла HCH (cch).
    Будем изменять его значения от 109.2 до 113.2 с шагом 0.2. В результате получим файл angle со значениями угла и соответствующими им энергиями. Зависимость энергии от угла была аппроксимирована параболой со значениями коэффициентов koeff2.

    График зависимости:





  7. Проделаем аналогичные операции для торсионного угла d3.
    Будем изменять его значения от -180 до 180 с шагом 12. В результате получим файл tors со значениями угла и соответствующими им энергиями. Зависимость энергии от угла:





    Количество минимумов: 3 (правый и левый совпадают).