Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.fbb.msu.ru/FBB/year_08/term6/task_4.html
Дата изменения: Wed Mar 16 17:13:14 2011
Дата индексирования: Tue Oct 2 07:43:34 2012
Кодировка: UTF-8
Task 4, term 6

Занятие 4.

Отчет по заданию должен появиться на сайте к следующему занятию.
Необходимые сведения о работе с GAMESS см. здесь.
Сведения о работе с Gnuplot см. здесь.
Введения о скриптовании в Bash здесь.
Ведения о awk здесь.
Вся работа по расчетам будет проходить на kodomo через терминал putty, а для работы с графическим выводом Gnuplot понадобится Xming.
Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчетов. Пример подобной работы представлен в статье о разработке поля gaff.
  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. Bash это интерпретатор командной строки который автоматически запускается при присоединении к kodomo.
  2. Составьте файл-заготовку для размножения, пусть имя файла будет et.inp. Для этого к координатам добавьте шапку для dft из предыдущего практикума. Только надо изменить информацию о типе входных координат: замените COORD=CART на COORD=ZMT. Проверьте работает ли Ваш файл-заготовка, т.е. запустите GAMESS если выходной файл не содержит ошибок, то переходим к следующему пункту.

  3. Теперь давайте создадим текстовый файл скрипта make_b.bash со следующим содержанием:

    #!/bin/bash
    ### делаем цикл от -10 до 10 ##### 
    
    for i in {-10..10}; do 
    
    #### нам надо рассчитать новую длину связи #####
    #### с шагом 0.02 ангстрема,               #####
    #### воспользуемся калькулятором bc        #####
    #### и результат поместим в переменную nb  #####
    
         nb=$(echo "scale=5; 1.001 + $i/50" | bc -l)
    
    #### пролистаем файл et.inp и заменим указание переменной ###
    #### на новое значение и пере направим результат в файл   ###
    
         sed "s/cc=1.001/cc=$nb/" et.inp  > b_${i}.inp
    
    done
    
    Поправьте скрипт так, что бы стартовая длина изменяемой связи соответствовала файлу et.inp. Запустите Ваш скрипт :
     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
  5. У вас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel. Вам предлагается сделать это в gnuplot.
    Запустите Xming->XLaunch. Выберите тип расположения окон, удобно использовать Multiple windows. Next. Выберите Start Program. Run Remote-> Putty. Дальше все как обычно: kodomo, username.
    Перейдите в рабочую директорию. Запустите Gnuplot:

    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 
    Сохраните значения коэффициентов. Постройте графики функции и значений энергии из Gamess.
    plot "bond", f(x)
    Сохраните изображение для отчета. В отчете обсудите причину не точного совпадения точек и функции.
  6. Проделайте аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2. Сохраните полученные коэффициенты в отчет. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчета. Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.

  7. Проделайте аналогичные операции для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12. Подгонку проводить не надо. Сохраните изображение и укажите в отчете количество минимумов функции.

  8. * Увеличьте шаг до 0.1 ангстрема при расчете связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.