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

Оптимизированная структура этана в виде z-матрицы c добавленной шапкой для DFT (тип входных координат - ZMT):

 $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
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

Проверим корректность файла, запустив GAMESS:

In [2]: 
%%bash
gms et.inp  1 >& et.log

Создадим на основе et.inp 20 файлов с разными длинами связи cc (шаг 0.02А):

In [9]: 
%%bash
for i in {-10..10}
do
    nb=$(echo "scale=5; 1.001 + $i*0.02" | bc -l)
    sed "s/cc=1.52986/cc=$nb/" et.inp  > b_${i}.inp
done

Запустим рассчет энергии с помощью GAMESS для получившихся файлов:

In []: 
%%bash
for i in {-10..10}; do gms b_${i}.inp 1 >& b_${i}.log; done

Запишем в один файл получившиеся значения энергий:

In [14]: 
%%bash
for i in {-10..10}
do 
    nb=$(echo "scale=5; 1.001 + $i*0.02" | bc -l)
    echo -n "$nb    "
    awk '/TOTAL ENERGY =/{print $4}'  b_${i}.log
done > bond
In [30]: 
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [33]: 
xx=[];yy=[]
for line in open("bond",'r'):
    x,y=map(float,line.split())
    xx.append(x);yy.append(y)

plt.scatter(xx, yy)
plt.xlabel("C-C length")
plt.ylabel("Energy")
Out[33]:
<matplotlib.text.Text at 0x110c9bb90>

Попробуем апроксимировать зависимость функцией \(f(x)=a+k(x-b)^2\):

In [37]: 
fit=np.polyfit(xx,yy,2)
xfit=[ min(xx)+(max(xx)-min(xx))/100*i for i in range(100) ]

plt.scatter(xx, yy)
plt.plot(xfit,np.polyval(fit, xfit),'r')
Out[37]:
[<matplotlib.lines.Line2D at 0x1163d9a10>]
In [38]: 
print fit

[  7.24461068 -17.38404908 -69.20171358]

\(f(x)=7.245x^2-17.384x-69.202\)

\(k=7.245\)

\(b=17.384 / 2k=1.200\)

\(a=-69.202-kb^2=-79.635\)

Проделаем то же самое для валентного угла cchv (значения изменяются от 109.2 до 113.2):

In [45]: 
%%bash
for i in {-10..10}
do
    nb=$(echo "scale=5; 111.2 + $i*0.2" | bc -l)
    sed "s/cchv=111.200/cchv=$nb/" et.inp  > v_${i}.inp
done
In [48]: 
%%bash
for i in {-10..10}; do gms v_${i}.inp 1 >& v_${i}.log; done
In [63]: 
%%bash
for i in {-10..10}
do 
    nb=$(echo "scale=5; 111.2 + $i*0.2" | bc -l)
    echo -n "$nb    "
    awk '/TOTAL ENERGY =/{print $4}'  v_${i}.log
done > vangle
In [62]: 
xx=[];yy=[]
for line in open("vangle",'r'):
    x,y=map(float,line.split())
    xx.append(x);yy.append(y)

plt.clf()
plt.plot(xx, yy,".")
plt.xlabel("HCC angle")
plt.ylabel("Energy")

fit=np.polyfit(xx,yy,2)
xfit=[ min(xx)+(max(xx)-min(xx))/100*i for i in range(100) ]

plt.plot(xfit,np.polyval(fit, xfit),'r')
Out[62]:
[<matplotlib.lines.Line2D at 0x116ecea50>]

Проделаем то же самое для торсионного угла d3 (значения изменяются от -180 до 180 c шагом 12):

In [67]: 
%%bash
for i in {-15..15}
do
    nb=$(echo "scale=5; 0 + $i*12" | bc -l)
    sed "s/d3=180/d3=$nb/" et.inp  > t_${i}.inp
done
In [48]: 
%%bash
for i in {-15..15}; do gms t_${i}.inp 1 >& t_${i}.log; done
In [70]: 
%%bash
for i in {-15..15}
do 
    nb=$(echo "scale=5; 0 + $i*12" | bc -l)
    echo -n "$nb    "
    awk '/TOTAL ENERGY =/{print $4}'  t_${i}.log
done > torsion
In [71]: 
xx=[];yy=[]
for line in open("torsion",'r'):
    x,y=map(float,line.split())
    xx.append(x);yy.append(y)

plt.clf()
plt.plot(xx, yy,".")
plt.xlabel("Torsion angle")
plt.ylabel("Energy")
Out[71]:
<matplotlib.text.Text at 0x116d6da10>

Мы наблюдаем три минимума (которые соответствуют, по-видимому, заторможенным конформациям).