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

Задание №1

Было решено автоматизировать процесс получения соответствующего числа файла через скрипт bash. Для начала был создан файл et.inp со следующим заголовком: 

 $CONTRL COORD=ZMT UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY $END
 $BASIS  GBASIS=N31 NGAUSS=6
  POLAR=POPN31 NDFUNC=1 $END
 $GUESS  GUESS=HUCKEL $END
 $system mwords=2 $end 
 $DATA

Если запустить GAMESS с таким файлом, то в выходном ошибок не будет, а это значит, что все ок! Далее был сделан bash файл, как описано в задании. После его запуска получился 21 файл (надо проследить, что файл находится в кодировке UNIX, иначе будут ругаться на концы строк). В итоге был получен файл bond в котором хранится следующая информация: 

1.32986    -79.7339315289
1.34986    -79.7406173737
1.36986    -79.7462868578
1.38986    -79.7510331813
1.40986    -79.7549412921
1.42986    -79.7580886263
1.44986    -79.7605457819
1.46986    -79.7623771323
1.48986    -79.7636413832
1.50986    -79.7643920816
1.52986    -79.7646780790
1.54986    -79.7645439543
1.56986    -79.7640303992
1.58986    -79.7631745699
1.60986    -79.7620104053
1.62986    -79.7605689179
1.64986    -79.7588784564
1.66986    -79.7569649422
1.68986    -79.7548520851
1.70986    -79.7525615748
1.72986    -79.7501132582
А дальше дело техники! Построенная зависимость энергии от длины C-C связи молекулы этана:
Параметры f(x): 
 
a = -79.7652 
k = 0.563608 
b = 1.55432
Как видно, функция далеко не совсем точно совпадает с полученными точками. Почему так? Ну, возможно, реальная зависимость энергии связи от ее длины описывается куда более сложной функцией!

Задание №2

Далее нужно было проделать аналогичные процедуры, но уже для валентного угла HCH. Вот соответствующие значения торсионного угла и энергии: 

109.200    -79.7645100964
109.400    -79.7645396600
109.600    -79.7645663953
109.800    -79.7645902988
110.000    -79.7646113672
110.200    -79.7646295973
110.400    -79.7646449862
110.600    -79.7646575313
110.800    -79.7646672302
111.000    -79.7646740809
111.200    -79.7646780815
111.400    -79.7646792309
111.600    -79.7646775278
111.800    -79.7646729717
112.000    -79.7646655623
112.200    -79.7646552997
112.400    -79.7646421847
112.600    -79.7646262180
112.800    -79.7646074011
113.000    -79.7645857358
113.200    -79.7645612241
Теперь непосредственно график.
Параметры f(x): 
 
a = -79.7647 
k = 3.56076e-05 
b = 111.38
Тут можно сказать, что апроксимация выполнена хорошо.

Задание №3

Теперь что касается торсионного угла. Проводим аналогичные процедуры, и получаем соответствующие значения для угла и энергии: 

-180    -79.7646780790
-168    -79.7642336928
-156    -79.7630628377
-144    -79.7616148156
-132    -79.7604407798
-120    -79.7599827599
-108    -79.7604407798
-96    -79.7616148156
-84    -79.7630628377
-72    -79.7642336928
-60    -79.7646780790
-48    -79.7642336928
-36    -79.7630628377
-24    -79.7616148156
-12    -79.7604407798
0    -79.7599827599
12    -79.7604407798
24    -79.7616148156
36    -79.7630628377
48    -79.7642336928
60    -79.7646780790
72    -79.7642336928
84    -79.7630628377
96    -79.7616148156
108    -79.7604407798
120    -79.7599827599
132    -79.7604407798
144    -79.7616148156
156    -79.7630628377
168    -79.7642336928
180    -79.7646780790
Вот график без подгонки.
Здесь видно 3 минимума, а не 4 (ибо значения для -180 и +180 совпадают - надо быть внимательнее!).