Домашняя работа 4. Вычисление параметров молекулярной механики
Основная рабочая директория с файлами: H:83
Создаю входной файл для GAMESS "et.inp" и запускаю GAMESS:
In [2]:
with open('et.inp','w') as f:
f.write(""" $CONTRL COORD=ZMT UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE $END
$BASIS GBASIS=N31 NGAUSS=6 $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""")
In [3]:
%%bash
gms et.inp 1 >& et.log
GAMESS работает с файлом, можно благополучно заменять переменные из конца файла.
Изменение длины связи cc
In [4]:
with open("make_b.sh",'w') as f:
f.write("""#!/bin/bash
cd bond
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""")
In [5]:
%%bash
bash ./make_b.sh
In [6]:
with open("run_gamess.sh",'w') as f:
f.write("""#!/bin/bash
cd bond
for i in {-10..10}; do
gms b_${i}.inp 1 > b_${i}.log
done""")
In [7]:
%%bash
bash ./run_gamess.sh
In [8]:
with open("parse_gamess_logs.sh",'w') as f:
f.write("""#!/bin/bash
cd bond
echo "Index\tcc\tEnergy1\t\tEnergy2"
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i/50" | bc -l)
echo -n "$i\t"
echo -n "$nb\t"
awk '/TOTAL ENERGY =/{printf "%s\t",$4}' b_${i}.log
echo ""
done""")
In [11]:
%%bash
bash ./parse_gamess_logs.sh > bond_test.txt
bash ./parse_gamess_logs.sh
Пробую построить и аппроксимировать полученные энергии. (Ссылки на хорошие гайды ipython notebook: по matplotlib и по fit-моделям в ipython notebook.
In [12]:
%matplotlib inline
In [13]:
import matplotlib.pyplot as plt
from pylab import *
from scipy.optimize import curve_fit
In [14]:
with open("bond_test.txt",'r') as f:
inf = f.readlines()[1:]
full_l = [[float(y) for y in x.split()] for x in inf]
ind_l = [int(x.split()[0]) for x in inf]
cc = [float(x.split()[1]) for x in inf]
e1 = [float(x.split()[2]) for x in inf]
e2 = [float(x.split()[3]) if (len(x.split())>3) else float(x.split()[2]) for x in inf ]
In [15]:
figure()
plot(cc, e1, 'k.')
xlabel('cc')
ylabel('energy1')
show()
In [16]:
figure()
plot(cc, e2, 'k.')
xlabel('cc')
ylabel('energy2')
show()
Видимо, значение второй энергии из файла - не нужно нам, т.к. представляет собой какие-то случайные значения.
In [17]:
def fitfunc(x, a, b, k):
return a + k*np.square(x-b)
#a+k(x-b)^2
fitpars, covmat = curve_fit(fitfunc,cc,e1)
variances = covmat.diagonal()
std_devs = np.sqrt(variances)
figure()
plot(cc, e1, 'k.')
plot(cc, fitfunc(cc,fitpars[0], fitpars[1],fitpars[2]), 'r-')
xlabel('cc')
ylabel('energy1')
show()
In [18]:
print fitpars
Видим неточное совпадение точек и функции.
Изменение валентного угла CCH
Диапазон изменения: от 109.2 до 113.2.
In [19]:
with open("make_angle.sh",'w') as f:
f.write("""#!/bin/bash
cd angle
for i in {0..20}; do
nb=$(echo "scale=5; 106.2 + $i*0.5" | bc -l)
sed "s/cch=111.200/cch=$nb/" ../et.inp > a_${i}.inp
done""")
with open("run_gamess_angle.sh",'w') as f:
f.write("""#!/bin/bash
cd angle
for i in {0..20}; do
gms a_${i}.inp 1 > a_${i}.log
done""")
with open("parse_gamess_logs_angle.sh",'w') as f:
f.write("""#!/bin/bash
cd angle
echo "Index\tcch\tEnergy1\t\tEnergy2"
for i in {0..20}; do
nb=$(echo "scale=5; 106.2 + $i*0.5" | bc -l)
echo -n "$i\t"
echo -n "$nb\t"
awk '/TOTAL ENERGY =/{printf "%s\t",$4}' a_${i}.log
echo ""
done""")
In [20]:
%%bash
bash ./make_angle.sh
bash ./run_gamess_angle.sh
bash ./parse_gamess_logs_angle.sh > angle_test.txt
bash ./parse_gamess_logs_angle.sh
In [21]:
with open("angle_test.txt",'r') as f:
inf = f.readlines()[1:]
full_l = [[float(y) for y in x.split()] for x in inf]
ind_l = [int(x.split()[0]) for x in inf]
cch = [float(x.split()[1]) for x in inf]
e1_a = [float(x.split()[2]) for x in inf]
e2_a = [float(x.split()[3]) if (len(x.split())>3) else float(x.split()[2]) for x in inf ]
In [22]:
figure()
plot(cch, e1_a, 'k.')
xlabel('cch')
ylabel('energy1')
show()
In [23]:
figure()
plot(cch, e2_a, 'k.')
xlabel('cch')
ylabel('energy2')
show()
In [25]:
def fitfunc(x, a, b, k):
return a + k*np.square(x-b)
#a+k(x-b)^2
fitpars_a, covmat = curve_fit(fitfunc,cch,e1_a)
figure()
plot(cch, e1_a, 'k.')
plot(cch, fitfunc(cch,fitpars_a[0], fitpars_a[1],fitpars_a[2]), 'r-')
xlabel('cch')
ylabel('energy1')
show()
In [26]:
print fitpars_a
Изменение торсионного угла d3
Диапазон изменения: от -180 до 180 c шагом 12.
In [27]:
with open("make_torsion.sh",'w') as f:
f.write("""#!/bin/bash
cd torsion
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""")
with open("run_gamess_torsion.sh",'w') as f:
f.write("""#!/bin/bash
cd torsion
for i in {-15..15}; do
gms t_${i}.inp 1 > t_${i}.log
done""")
with open("parse_gamess_logs_torsion.sh",'w') as f:
f.write("""#!/bin/bash
cd torsion
echo "Index\td3\tEnergy1\t\tEnergy2"
for i in {-15..15}; do
nb=$(echo "scale=5; 0 + $i*180/15" | bc -l)
echo -n "$i\t"
echo -n "$nb\t"
awk '/TOTAL ENERGY =/{printf "%s\t",$4}' t_${i}.log
echo ""
done""")
In [28]:
%%bash
bash ./make_torsion.sh
bash ./run_gamess_torsion.sh
bash ./parse_gamess_logs_torsion.sh > torsion_test.txt
bash ./parse_gamess_logs_torsion.sh
In [29]:
with open("torsion_test.txt",'r') as f:
inf = f.readlines()[1:]
full_l = [[float(y) for y in x.split()] for x in inf]
ind_l = [int(x.split()[0]) for x in inf]
d3 = [float(x.split()[1]) for x in inf]
e1_t = [float(x.split()[2]) for x in inf]
e2_t = [float(x.split()[3]) if (len(x.split())>3) else float(x.split()[2]) for x in inf ]
In [30]:
figure()
plot(d3, e1_t, 'k.')
xlabel('d3')
ylabel('energy1')
show()
In [31]:
figure()
plot(d3, e2_t, 'k.')
xlabel('d3')
ylabel('energy2')
show()
Вижу 4 минимума и 3 максимума функции энергии в зависимости от торсионного угла. energy2 из файла я не рассматриваю.
energy1 представляет собой колебательную функцию, как и можно ожидать для изменения торсионного угла. Максимумы энергии - заслоненные конформации, минимумы - заторможенные.
Увеличение шага длины связи
Итак, пробую увеличить шаг до 0.1 (было 0.02):
In [32]:
with open("make_b1.sh",'w') as f:
f.write("""#!/bin/bash
cd bond1
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i*0.1" | bc -l)
sed "s/cc=1.52986/cc=$nb/" ../et.inp > b1_${i}.inp
done""")
with open("run_gamess_bond1.sh",'w') as f:
f.write("""#!/bin/bash
cd bond1
for i in {-10..10}; do
gms b1_${i}.inp 1 > b1_${i}.log
done""")
with open("parse_gamess_logs_bond1.sh",'w') as f:
f.write("""#!/bin/bash
cd bond1
echo "Index\tcch\tEnergy1\t\tEnergy2"
for i in {-10..10}; do
nb=$(echo "scale=5; 1.52986 + $i*0.1" | bc -l)
echo -n "$i\t"
echo -n "$nb\t"
awk '/TOTAL ENERGY =/{printf "%s\t",$4}' b1_${i}.log
echo ""
done""")
In [33]:
%%bash
bash ./make_b1.sh
bash ./run_gamess_bond1.sh
bash ./parse_gamess_logs_bond1.sh > bond1_test.txt
bash ./parse_gamess_logs_bond1.sh
In [34]:
with open("bond1_test.txt",'r') as f:
inf = f.readlines()[1:]
full_l = [[float(y) for y in x.split()] for x in inf]
ind_l = [int(x.split()[0]) for x in inf]
cc1 = [float(x.split()[1]) for x in inf]
e1_b1 = [float(x.split()[2]) for x in inf]
e2_b1 = [float(x.split()[3]) if (len(x.split())>3) else float(x.split()[2]) for x in inf ]
In [35]:
figure()
plot(cc1, e1_b1, 'k.')
xlabel('cc')
ylabel('energy1')
show()
In [36]:
figure()
plot(cc1[1:], e2_b1[1:], 'k.')
xlabel('cc')
ylabel('energy2')
show()
Из графика energy1 видно, что вид зависимости похож на гиперболу (по крайней мере в первом приближении).
In [1]:
%%bash
# Просто часто используемая функция
rm /home/students/y11/agalicina/gamess-scratch/*