Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.cmm.msu.ru/~lynx/term6/text/molmech.html
Дата изменения: Wed Mar 14 02:54:12 2012
Дата индексирования: Tue Oct 2 17:11:31 2012
Кодировка: Windows-1251
MolMech

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

    Определение зависимости энергии молекулы этана от длины связи C-C

  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 с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не ее значением.
  2. Автоматизируем процесс с помощью скрипта на bash. Bash это интерпретатор командной строки, который автоматически запускается при присоединении к kodomo. Составим файл-заготовку для размножения, пусть имя файла будет et.inp. Для этого к координатам добавим шапку для dft из предыдущего практикума. Только надо изменить информацию о типе входных координат: заменить COORD=CART на COORD=ZMT.
  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 стaвим запуск Gamess:
    gms b_${i}.inp 1 > b_${i}.log
    Для извлечения значения энергий из log файла закомментируем в скрипте запуск Gamess. Затем добавим вызов awk:
    echo -n "$nb " (pаспечатаем переменную nb и несколько пробелов без переноса строки);
    awk '/TOTAL ENERGY =/{print $4}' b_${i}.log (ищем строчку с TOTAL ENERGY и печатаем 4ое поле считая что поля разделены пробелами).
    После запуска скрипта перенаправим в файл две колонки цифр (зависимость энергии от длины связи):
    bash ./make_b.bash > bond
  5. Зависимость энергии молекулы от длины одной связи можно в 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
    Сохраним значения коэффициентов:
    
    
    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%)
    
    
    correlation matrix of the fit parameters:
    
                   a      k      b
    a               1.000
    k              -0.724  1.000
    b               0.171 -0.413  1.000
    

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



    Cамая большая ошибка для k. Возможно, из-за этого полученная функция и точки совпадают не абсолютно.

    Определение зависимости энергии молекулы этана от значения валентного угла HCH

  6. Для валентного угла HCH были проведены аналогичные операции, причем его занчения изменялись от 109.2 до 113.2 (файл inp,скрипт, значения углов HCH).
    Коэффициенты f(x), подогнанные gnuplot (сама функция и стартовые значения коэффициентов аналогично предыдущему заданию):
    
    
    Final set of parameters            Asymptotic Standard Error
    =======================            ==========================
    
    a               = -79.7647         +/- 1.21e-08     (1.517e-08%)
    k               = 3.56076e-05      +/- 6.229e-09    (0.01749%)
    b               = 111.38           +/- 9.954e-05    (8.937e-05%)




  7. Определение зависимости энергии молекулы этана от значения торсионного угла d3

  8. Для торсионного угла в3 были проведены аналогичные операции, причем его занчения изменялись от 180 до -180.
    По полученным данным в gnuplot был построен график. Он похож на синусоиду, поэтому эта зависимость была аппроксимирована следующей функцией:
    f(x)=k*cos(n*x)-a
    где k - константа жесткости, a - смещение, n - коэффициент растяжения/сжатия. Их начальные значения:
    a = -80
    n = 3
    k = 1
    После подгонки коэффициенты приобрели следующий вид:
    
    
    Final set of parameters            Asymptotic Standard Error
    =======================            ==========================
    
    a               = -79.7623         +/- 6.577e-07    (8.245e-07%)
    n               = 3.00014          +/- 0.0002247    (0.007491%)
    k               = 0.00234519       +/- 8.891e-07    (0.03791%)
    
    

    График, полученный в gnuplot:



    Мы наблюдаем 3 минимума, так как крайние точки совпадают (угол -180 равне углу 180).
  9. Определение зависимости энергии молекулы этана от длины связи C-C при увеличении шага до 0.1 ангстрема

  10. скрипт
    полученые значения
    Значения коэффициентов:
    
    Final set of parameters            Asymptotic Standard Error
    =======================            ==========================
    
    a               = -79.7647         +/- 5.444e-05    (6.825e-05%)
    k               = 0.53427          +/- 0.01117      (2.091%)
    b               = 1.53859          +/- 0.0005931    (0.03855%)
    
    
    correlation matrix of the fit parameters:
    
                   a      k      b      
    a               1.000 
    k              -0.734  1.000 
    b               0.126 -0.308  1.000 




    Аппроксимация получилась ровно такая же (и функция, и значения параметров), как и в самом первом задании. Но связь получилась короче и энергия повыше (в первом задании она была в диапазоне [-79,77...-79,73], здесь же [-79,765...-79,758]).

© Anastasia Maslova, 2012