Практикум №4. Вычисление ковалентных и нековалентных параметров для молекулярной механики¶
1-5. Константы для СС-связи этана¶
Оптимизированная структура этана в виде z-матрицы, подготовленная для расчета энергии по теории функционала плотности:
$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
ethane
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
Создание файлов, в которых длина СС-связи этана различается на 0.02 А, и расчет энергии этана по теории функционала плотности с помощью GAMESS:
%%bash
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)
sed "s/cc=1.52986/cc=$nb/" ethane_zmt_dft.inp > cc_bond_${i}.inp
gms cc_bond_${i}.inp 1 > cc_bond_${i}.log
done
Построение графика зависимости энергии этана от длины СС-связи, аппроксимация полученной зависимости функцией \(f(x)=a+k(x-b)^2\) и сохранение координат в файл:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
cc_bonds = list()
energies_1 = list()
cc_bond = open('cc_bond.txt', 'w')
for i in range(-10, 11):
log = open('cc_bond_%s.log'%i, 'r')
for line in log:
if line.strip().startswith('INPUT CARD>cc='):
bond = line.split('=')[1].strip()
cc_bonds.append(float(bond))
cc_bond.write(bond+' ')
if line.strip().startswith('TOTAL ENERGY'):
energy = line.split('=')[1].strip()
energies_1.append(float(energy))
cc_bond.write(energy+'\n')
log.close()
cc_bond.close()
[a0, b0, c0] = np.polyfit(cc_bonds, energies_1, 2)
print 'a0 = %s, b0 = %s, c0 = %s'%(a0, b0, c0)
# a0 = k
# b0 = -2kb
# c0 = a+kb^2
k = a0
b = b0/(-2*a0)
a = c0-a0*(b0/(-2*a0))**2
print 'k = %s, b = %s, a = %s'%(k, b, a)
new_energies_1 = list()
for x in cc_bonds:
y = a+k*(x-b)**2
new_energies_1.append(y)
plt.plot(cc_bonds, new_energies_1, label='f(x)')
plt.plot(cc_bonds, energies_1, color='lightgreen', marker='o', ls='*', label='CC-bond')
plt.legend(loc='upper center')
plt.xlabel("CC-bond length (Angstrom)")
plt.ylabel("Ethane energy (Hartree)")
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
from IPython.display import FileLink
FileLink('cc_bond.txt')
На графике выше видно, что зависимость энергии этана от длины СС-связи не является квадратичной.
6. Константы для НСН-угла этана¶
Создание файлов, в которых величина HСH-угла этана различается на 0.2°, и расчет энергии этана по теории функционала плотности с помощью GAMESS:
%%bash
for i in {-10..10}; do
nb=$(echo "scale=1; 111.2 + $i/5" | bc -l)
sed "s/111.200/$nb/" ethane_zmt_dft.inp > hch_angle_${i}.inp
gms hch_angle_${i}.inp 1 > hch_angle_${i}.log
done
Построение графика зависимости энергии этана от величины НСН-угла, аппроксимация полученной зависимости функцией \(f(x)=a+k(x-b)^2\) и сохранение координат в файл:
hch_angles = list()
energies_2 = list()
hch_angle = open('hch_angle.txt', 'w')
for i in range(-10, 11):
log = open('hch_angle_%s.log'%i, 'r')
for line in log:
if line.strip().startswith('INPUT CARD>cch='):
angle = line.split('=')[1].strip()
hch_angles.append(float(angle))
hch_angle.write(angle+' ')
if line.strip().startswith('TOTAL ENERGY'):
energy = line.split('=')[1].strip()
energies_2.append(float(energy))
hch_angle.write(energy+'\n')
log.close()
hch_angle.close()
[a0, b0, c0] = np.polyfit(hch_angles, energies_2, 2)
print 'a0 = %s, b0 = %s, c0 = %s'%(a0, b0, c0)
k = a0
b = b0/(-2*a0)
a = c0-a0*(b0/(-2*a0))**2
print 'k = %s, b = %s, a = %s'%(k, b, a)
new_energies_2 = list()
for x in hch_angles:
y = a+k*(x-b)**2
new_energies_2.append(y)
plt.plot(hch_angles, new_energies_2, label='f(x)')
plt.plot(hch_angles, energies_2, color='lightgreen', marker='o', ls='*', label='HCH-angle')
plt.legend(loc='upper center')
plt.xlabel("HCH-angle value (Degree)")
plt.ylabel("Ethane energy (Hartree)")
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
FileLink('hch_angle.txt')
На графике выше видно, что зависимость энергии этана от величины НСН-угла хорошо аппроксимируется квадратичной функцией.
7. Торсионный d3-угол этана¶
Создание файлов, в которых величина d3-угла этана различается на 12°, и расчет энергии этана по теории функционала плотности с помощью GAMESS:
%%bash
for i in {-30..0}; do
nb=$(echo "scale=0; 180 + $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" ethane_zmt_dft.inp > d3_angle_${i}.inp
gms d3_angle_${i}.inp 1 > d3_angle_${i}.log
done
Построение графика зависимости энергии этана от величины d3-угла и сохранение координат в файл:
angles_d3 = list()
energies_3 = list()
d3_angle = open('d3_angle.txt', 'w')
for i in range(-30, 1):
log = open('d3_angle_%s.log'%i, 'r')
for line in log:
if line.strip().startswith('INPUT CARD>d3='):
angle = line.split('=')[1].strip()
angles_d3.append(float(angle))
d3_angle.write(angle+' ')
if line.strip().startswith('TOTAL ENERGY'):
energy = line.split('=')[1].strip()
energies_3.append(float(energy))
d3_angle.write(energy+'\n')
log.close()
d3_angle.close()
plt.plot(angles_d3, energies_3, color='lightgreen', marker='o', ls='*', label='d3-angle')
plt.legend(loc='upper center')
plt.xlabel("d3-angle value (Degree)")
plt.ylabel("Ethane energy (Hartree)")
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
FileLink('d3_angle.txt')
На графике выше видно, что минимуму энергии этана соответствуют три торсионных d3-угла: -60°, 60° и 180°(-180°).