Вычисление параметров для молекулярной механики
import sys
header=""" $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"""
with open('et.inp', 'w') as et:
et.write(header)
Полученный файл запустили в GAMESS, для проверки отсутствия ошибок в выходном файле.
%%bash
gms et.inp > et.log
Ошибок не было получено,следовательно, нужно создать с помощью скрипта bash 21 файла с различными значениями переменной cc.
%%bash
bash ./make_b.bash
Модернизируем скрипт для запуска расчетов для всех этих файлов, и запустим его.
%%bash
bash ./game_make.bash
Затем изменим примененный ранее скрипт так, чтобы на экран выводились значения энергии и длины связей, полученных при расчетах GAMESS.
%%bash
bash ./n_make.bash > bond.txt
bash ./n_make.bash
Визуализировали полученную зависимость с помощью matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
a=np.loadtxt('bond.txt')
x=a[:,0]
y=a[:,1]
plt.plot(x, y,'k.')
plt.xlabel("C-C bond length")
plt.ylabel("Total energy")
Теперь проведем описание полученной зависимости функцией:
a=np.loadtxt('bond.txt')
x=a[:,0]
y=a[:,1]
coefficients = np.polyfit(x, y, 2)
polynomial = np.poly1d(coefficients)
xs = np.linspace(min(x),max(x),50)
ys = polynomial(xs)
plt.plot(x, y, 'k.')
plt.plot(xs, ys)
plt.xlabel("C-C bond length")
plt.ylabel("Total energy")
print polynomial
Наблюдаются неточности в описываемой функции.
Далее проведем аналогичные операции для валентного угла HCH, изменяя его значения от 109.2 до 113.2. Для этого напишем аналогичные предыдущим скрипты.
with open("make_a.sh",'w') as file:
file.write("""
for i in {-10..10}; do
nb=$(echo "scale=5; 111.200 + $i/5" | bc -l)
sed "s/cchv=111.200/cchv=$nb/" et.inp > a_${i}.inp
done""")
with open("game_a.sh",'w') as file:
file.write("""
for i in {-10..10}; do
nb=$(echo "scale=5; 111.200 + $i/5" | bc -l)
sed "s/cchv=111.200/cchv=$nb/" et.inp > a_${i}.inp
gms a_${i}.inp 1>& a_${i}.log
done""")
with open("n_make_a.sh",'w') as file:
file.write("""
for i in {-10..10}; do
nb=$(echo "scale=5; 111.200 + $i/5" | bc -l)
sed "s/cchv=111.200/cchv=$nb/" et.inp > a_${i}.inp
echo -n "$nb "
awk '/TOTAL ENERGY =/{print $4}' a_${i}.log
done""")
%%bash
bash ./make_a.sh
%%bash
bash ./game_a.sh
%%bash
bash ./n_make_a.sh > angle.txt
less angle.txt
a=np.loadtxt('angle.txt')
x=a[:,0]
y=a[:,1]
coefficients = np.polyfit(x, y, 2)
polynomial = np.poly1d(coefficients)
xs = np.linspace(min(x),max(x),50)
ys = polynomial(xs)
plt.plot(x, y, 'k.')
plt.plot(xs, ys)
plt.xlabel("HCH angle")
plt.ylabel("Total energy")
print polynomial
Все тоже самое сделаем для торсионного угла d3, изменяя его значения от -180 до 180 c шагом 12. Для этого напишем аналогичные предыдущим скрипты.
with open("make_d.sh",'w') as file:
file.write("""
for i in {-15..15}; do
nb=$(echo "scale=0; 0 + $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" et.inp > c_${i}.inp
done""")
with open("game_d.sh",'w') as file:
file.write("""
for i in {-15..15}; do
nb=$(echo "scale=0; 0 + $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" et.inp > c_${i}.inp
gms c_${i}.inp 1>& c_${i}.log
done""")
with open("n_make_d.sh",'w') as file:
file.write("""
for i in {-15..15}; do
nb=$(echo "scale=0; 0 + $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" et.inp > c_${i}.inp
gms c_${i}.inp 1>& c_${i}.log
echo -n "$nb"
awk '/TOTAL ENERGY =/{print $4}' c_${i}.log
done""")
%%bash
bash ./make_d.sh
%%bash
bash ./game_d.sh
%%bash
bash ./n_make_d.sh > d.txt
less d.txt
a=np.loadtxt('d.txt')
x=a[:,0]
y=a[:,1]
plt.plot(x, y, 'k.')
plt.xlabel("d3 torsion angle")
plt.ylabel("Total energy")
На графике отчетливо видны 3 минимума функции(-180 и 180 прядставляют собой одну и ту же точку.)