Была составлена заготовка для размножения:

 $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 было проверено, что она не содержит в себе ошибок.

3

Потом с помощью скрипта было создано еще 20 подобных файлов, в которых длина связи C-C изменялась с шагом 0.02 ангстрема:

In [4]:

%%bash
for i in {-10..10}; do 
    nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)
    sed "s/cc=1.52986/cc=$nb/" et.inp  > b_${i}.inp
done

В итоге получен 21 файл с разным значением переменной сс.

4

В скрипт были добавлены строчки для запуска gamess и поиска значений энергии:

In [12]:

%%bash
echo 'Length     Energy'
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)
sed "s/cc=1.52986/cc=$nb/" et.inp  > b_${i}.inp
gms b_${i}.inp 1 > b_${i}.log
echo -n "$nb    "
awk '/TOTAL ENERGY =/{print $4}' b_${i}.log

done

Length     Energy
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

5

Был построен график полученной зависимости. Затем с помощью gnuplot была проведена его аппроксимация параболой \(f(x)= a+k(x-b)^2\). Полученные коэффициенты:

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

На рисунке ниже представлены данные, полученные при помощи расчетов (точки), и функция с подогнанными коэффициентами. Видно, что аппроксимация квадратичной функцией неточная, соответственно, связь энергии молекулы от длины связей подчиняется какой-то другой зависимости (а именно потенциалу Леннард-Джонса, что будет видно позднее).

In [9]:

Image('cc.png')

Out[9]:

6

Аналогичные расчеты и построения были проделаны для валентного угла HCH. Зависимость энергий от угла получилась следующая:

In [16]:

%%bash
echo 'Angle     Energy'
for i in {-10..10}; do
nb=$(echo "scale=5; 111.2 + $i/5" | bc -l)
sed "s/cchv=111.2/cchv=$nb/" et.inp  > angle_${i}.inp
gms angle_${i}.inp 1 > angle_${i}.log
echo -n "$nb    "
awk '/TOTAL ENERGY =/{print $4}' angle_${i}.log

done

Angle     Energy
109.20000    -79.7645100655
109.40000    -79.7645396319
109.60000    -79.7645663700
109.80000    -79.7645902763
110.00000    -79.7646113475
110.20000    -79.7646295805
110.40000    -79.7646449722
110.60000    -79.7646575202
110.80000    -79.7646672219
111.00000    -79.7646740754
111.20000    -79.7646780790
111.40000    -79.7646792311
111.60000    -79.7646775309
111.80000    -79.7646729776
112.00000    -79.7646655711
112.20000    -79.7646553114
112.40000    -79.7646421992
112.60000    -79.7646262354
112.80000    -79.7646074214
113.00000    -79.7645857589
113.20000    -79.7645612501

In [15]:

Image('hch.png')

Out[15]:

7

Аналогичные операции (без аппроксимации) были проделаны для торсионного угла d3:

In [17]:

%%bash
echo 'Dihedral     Energy'
for i in {-15..15}; do
nb=$(echo "scale=5; $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" et.inp  > dih_${i}.inp
gms dih_${i}.inp 1 > dih_${i}.log
echo -n "$nb    "
awk '/TOTAL ENERGY =/{print $4}' dih_${i}.log

done

Dihedral     Energy
-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

In [18]:

Image('dih.png')

Out[18]:

Мы видим 3 минимума и 3 максимума функции (потому что -180 и +180 - один и тот же угол).

8

Зависимость энергии молекулы от длины связи С-С была расчитана еще раз, теперь с шагом 0.1 А. Результат:

In [19]:

%%bash
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i/10" | bc -l)
sed "s/cc=1.52986/cc=$nb/" et.inp  > ld_${i}.inp
gms ld_${i}.inp 1 > ld_${i}.log
echo -n "$nb    "
awk '/TOTAL ENERGY =/{print $4}' ld_${i}.log

done

.52986    -74.3435804461
.62986    -76.4967472763
.72986    -77.8283325186
.82986    -78.6344334655
.92986    -79.1198120262
1.02986    -79.4108122553
1.12986    -79.5828707009
1.22986    -79.6813083202
1.32986    -79.7339315289
1.42986    -79.7580886263
1.52986    -79.7646780790
1.62986    -79.7605689179
1.72986    -79.7501132582
1.82986    -79.7360890363
1.92986    -79.7202779104
2.02986    -79.7038189912
2.12986    -79.6874293573
2.22986    -79.6715494454
2.32986    -79.6564412757
2.42986    -79.6422534996
2.52986    -79.6290609314

В данном случае для аппроксимации надо использовать не квадратичную функцию, а потенциал Леннарда-Джонса. Для уравнения \(f(x) = 4e((\frac{k}{x-b})^{12} - \frac{k}{x-b}^6)) + a\) были подобраны следующие константы:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

e               = 0.224185         +/- 0.004382     (1.955%)
k               = 3.42195          +/- 0.007488     (0.2188%)
a               = -79.5393         +/- 0.003577     (0.004498%)
b               = -2.32613         +/- 0.009718     (0.4178%)

In [20]:

Image('ld.png')

Out[20]:

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

1

Сначала с помощью скриптов RED были определены точечные заряды молекулы. Для этого результат оптимизации был переформатирован в pdb, подготовлен с помощью Ante_RED.pl (мультиплетность молекулы = 1; заряд = 0, нас все устраивает) и наконец обработан с помощью RED-vIII.4.pl.

2

После этого был создан файл описания молекулы в формате пакета программ GROOMACS. В первых двух строчках задаются некоторые правила:

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes              0.5     0.8333

Затем задаются типы атомов и параметры для функции Леннарда-Джонса. Так как мы считали, что в случае этана Ван-дер-Ваальсово взаимодействие между атомами углерода разных молекул минимально, можно поставить для углерода некоторые параметры. Ван-дер-Ваальсов радиус известен, единственная переменная - epsilon водорода.

[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
H          1        1.008    0.0000  A   1.06908e-01  1.00000e-00
C          6        12.01    0.0000  A   3.39967e-01  3.59824e-01

Затем следует описание непосредственно молекулы. Указывается имя и то, что соседи через три не учитываются при расчете Ван-дер-Ваальсовых взаимодействий:

[ moleculetype ]
; Name            nrexcl
et            3

Затем добавляются атомы. Заряды были расчитаны в предыдущем задании.

[ atoms ]
;   nr  type  resnr  residue  atom   cgnr     charge       mass
     1   C      1    ETH      C1      1    -0.0189    12.01
     2   C      1    ETH      C2      2    -0.0155    12.01
     3   H      1    ETH      H1      3     0.0059    1.008
     4   H      1    ETH      H2      4     0.0059    1.008
     5   H      1    ETH      H3      5     0.0059    1.008
     6   H      1    ETH      H4      6     0.0059    1.008
     7   H      1    ETH      H5      7     0.0059    1.008
     8   H      1    ETH      H6      8     0.0059    1.008

Далее - описание связей.

[ bonds ]
;  ai    aj funct  b0       kb
     1   2   1  0.1554   150000.0 
     1   3   1  0.1085   180000.0
     1   4   1  0.1085   180000.0
     1   5   1  0.1085   180000.0
     2   6   1  0.1085   180000.0
     2   7   1  0.1085   180000.0
     2   8   1  0.1085   180000.0

Аналогично описываются углы и торсионные углы.

[ angles ]
;  ai    aj    ak funct  phi0   kphi
;around c1
    3     1     4     1   109.500    200.400
    3     1     5     1   109.500    200.400
    4     1     5     1   109.500    200.400
    3     1     2     1   109.500    400.400
    4     1     2     1   109.500    400.400
    5     1     2     1   109.500    400.400
;around c2
    1     2     6     1   109.500    400.400
    1     2     7     1   109.500    400.400
    1     2     8     1   109.500    400.400
    6     2     7     1   109.500    200.400
    6     2     8     1   109.500    200.400
    7     2     8     1   109.500    200.400

[ dihedrals ]
;  ai    aj    ak    al funct  t0           kt      mult
    3    1     2     6      1  0.0      0.62760     3
    3    1     2     7      1  0.0      0.62760     3
    3    1     2     8      1  0.0      0.62760     3
    4    1     2     6      1  0.0      0.62760     3
    4    1     2     7      1  0.0      0.62760     3
    4    1     2     8      1  0.0      0.62760     3
    5    1     2     6      1  0.0      0.62760     3
    5    1     2     7      1  0.0      0.62760     3
    5    1     2     8      1  0.0      0.62760     3

Затем - пары атомов, которые не должны считаться при расчете VdW.

[ pairs ]
;  ai    aj funct
   3  6
   3  7
   3  8
   4  6
   4  7
   4  8
   5  6
   5  7
   5  8

И, наконец, описание системы.

[ System ]
; eth
first one
[ molecules ]
;Name count
 et    38

et.top

3

Есть два состояния системы: первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем, второе имеет такую же плотность, как и жидкий этан. Надо провести короткое моделирование динамики каждой из этих систем, определить разницу в энергии VdW взаимодействий между ними и сравнить с энтальпией испарения этана (5.4 кДЖ/моль при T=25). При этом epsilon водорода нам неизвестна. Для этого создадим скрипт, который создает семь топологий с разными значениями epsilon, проводит для каждой системы молекулярную динамику с каждым файлом топологии и считает значения энергий:

In [26]:

%%bash
echo 'Epsilon    LJ(газ)   Coulomd(газ)   LJ(жидкость)  Coulomd(жидкость)'
for i in {1..7};do
    ep=$( echo "scale=5; 1/$i/$i/$i" | bc -l )
    sed "s/1.00000/$ep/" et.top > v_${i}.top
    grompp_d -f md -c box_big -p v_${i}.top -o vb_${i} -maxwarn 1 && mdrun_d -deffnm  vb_${i} -v
    grompp_d -f md -c box_38 -p v_${i}.top -o v_${i} -maxwarn 1 && mdrun_d -deffnm  v_${i} -v
    gas=`echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  vb_${i} -o eb_${i} | awk '(/LJ/ || /Coulomb/) {print $3}'`
    liq=`echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  v_${i} -o e_${i} | awk '(/LJ/ || /Coulomb/) {print $3}'`
        echo $ep    $gas    $liq > energy
    done
cat energy

Epsilon    LJ(газ)   Coulomd(газ)   LJ(жидкость)  Coulomd(жидкость)
1.00000   -0.000226809 3.89874e-05    -264.616 0.0017485
.12500    -0.00015805 3.89895e-05    -21.8114 0.000106461
.03703    -0.000141376 3.89889e-05    -5.01466 0.000213898
.01562    -0.000134447 3.89881e-05    -2.34411 0.000187274
.00800    -0.000130809 3.89874e-05    -1.93261 4.96359e-05
.00462    -0.000128615 3.89872e-05    -2.47454 5.15414e-05
.00291    -0.000127183 3.89875e-05    -3.3947 3.14326e-05

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

1

3

4

5

6

7

8

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

1

2

3