Вычисление параметров для молекулярной механики¶
In [23]:
from xmlrpclib import ServerProxy
import os, sys
from IPython.display import Image
os.environ['PATH']=os.environ['PATH']+'/home/preps/golovin/progs/bin'
os.environ['MOPAC_LICENSE'] ='/home/preps/golovin/progs/bin'
os.environ['LD_LIBRARY_PATH']='/home/preps/golovin/progs/bin'
In [3]:
%%bash
gms et.inp 1 >& et.log
Нам дана молекула этана. Cоздадим 21 разный файл для расчёта энергии в GAMESS с разными значениями по длине одной из связей, запустим расчёт для этих файлов.
In []:
%%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
gms b_${i}.inp 1 > b_${i}.log
done
Извлечём значение энергии из log файла.
In [21]:
%%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
#gms b_${i}.inp 1 > b_${i}.log
echo -n "$nb "
awk '/TOTAL ENERGY =/{print $4}' b_${i}.log
done > bond
Построим график зависимости энергии молекулы от длины одной связи с помощью Gnuplot
| 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 [24]:
Image(filename='cc.png')
Out[24]:
То же самое с matplot
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [26]:
a = np.loadtxt("bond")
x=a[:,0]
y=a[:,1]
print y
plt.scatter(x, y, alpha=0.5)
plt.show()
In [2]:
a = np.loadtxt("bond")
x_o=a[:,0]
y_o=a[:,1]
print "initial data:", y_o
#function is f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized params:", p1
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(1.3,1.8)
plt.show()
Для валентного угла HCH проварьируем значения от 109.2 до 113.2.
In []:
%%bash
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 > c_${i}.inp
gms c_${i}.inp 1 > c_${i}.log
echo -n "$nb "
awk '/TOTAL ENERGY =/{print $4}' c_${i}.log
done > cchv
In [17]:
a = np.loadtxt("cchv")
x_o=a[:,0]
y_o=a[:,1]
print "initial data:", y_o
#function is f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print "Optimized params:", p1
#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='magenta',alpha=0.5)
plt.xlim(109,113.4)
plt.show()
То же для торсионного угла d3. Рассмотрим все возможные занчения.
In []:
%%bash
for i in {-30..0}; do
nb=$(echo "scale=5; 180 + $i*12" | bc -l)
sed "s/d3=180/d3=$nb/" et.inp > d_${i}.inp
gms d_${i}.inp 1 > d_${i}.log
echo -n "$nb "
awk '/TOTAL ENERGY =/{print $4}' d_${i}.log
done > d3
In [5]:
a = np.loadtxt("d3")
x=a[:,0]
y=a[:,1]
print y
plt.scatter(x, y, alpha=0.5)
plt.xlim(-180,180)
plt.show()
Зависимость имеет три минимума, которые соответствуют незаслонённой конформации молекулы этана