Дана оптимизированная структура этана в виде 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.
Файл-заготовка для размножения: et.inp.
Для проверки файла запустим GAMESS.

Изменение длины сс-связи

Скрипт для размножения файла-заготовки: make_b1.bash. Изменяем длину связи сс от 1.32986 до 1.72986.
Запускаем скрипт:
bash ./make_b.bash
Мы получили 21 файл с разными значениями переменной сс.
Изменим скрипт так, чтобы запустить расчёт GAMESS: make_b2.bash.
Чтобы извлечь значение энергии из log файла, изменим скрипт ещё раз: make_b3.bash.
Запускаем скрипт и перенаправляем выход в файл bond:
bash ./make_b.bash > bond
В bond находятся значения энергии и соответствующие им длины связи.
В 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
Получились следующие значения коэффициентов:

a               = -79.7652         +/- 0.0004522    (0.000567%)
k               = 0.563608         +/- 0.02335      (4.142%)
b               = 1.55432          +/- 0.002455     (0.1579%)

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

Функция неточно совпадает с полученными точками, вероятно, потому, что реальная зависимость энергии связи от её длины описывается более сложной функцией.

Изменение валентного угла HCH

Аналогичные операции проделаем для валентного угла HCH, описываемого переменной cchv, изменяя его значения от 109.2 до 113.2.
Полученные в процессе файлы:
Скрипт на bash: make_b4.bash
Значения энергий: ecke
Аппроксимируем зависимость той же параболой:
f(x)=a + k*x*x - 2*k*x*b + k*b*b
Проведём подгонку коэффициентов и получим:

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%)

Аппроскимация прошла хорошо, все "экспериментальные" точки лежат на параболе.

Изменение торсионного угла d3

Аналогичные операции проделаем для торсионного угла d3, описываемого переменной d3, изменяя его значения от -180 до 180 c шагом 12.
Скрипт на bash: make_b5.bash
Значения энергий: tor
Аппроксимируем зависимость косинусом:
f(x)=a*cos(k*x*pi/180)+b
Проведём подгонку коэффициентов и получим:

a               = 0.00234519       +/- 8.891e-007   (0.03791%)
k               = 3.00014          +/- 0.0002247    (0.007491%)
b               = -79.7623         +/- 6.577e-007   (8.245e-007%)        

Количество минимумов - 3, так как значения d3 180 и -180 - это одна и та же точка.

---

©Andreeva_2010