Cеми-эмпирические вычисления: Mopac


Пункт 1. Модуль pybel

In [3]:
import pybel
In [4]:
mol=pybel.readstring('smi','c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3')
mol.addh()
mol.make3D(steps=10000)
In [6]:
mol.write(format='pdb',filename='pr5_1.pdb',overwrite=True)

Рис. 1. Структура порфирина, полученная модулем pybel.

Порфирин получается неправильный: 4 водорода у центральных атомов кислорода и изогнутый. Правильный: плоский и 2 атома водорода. Попробовали другим способом.

Пункт 2. MOPAC способ PM6/АМ1

In [11]:
mop=mol.write(format='mopin',opt={'k':'PM6 CHARGE=%d' % mol.charge})
In [12]:
mop=mol.write(format='mopin',filename='my.mop',opt={'k':'PM6 CHARGE=%d' % mol.charge},overwrite=True)

В my.mop заряд молекулы равен нулю, поэтому перейдем в способу АМ1

In [14]:
mop_am1=mol.write(format='mopin',opt={'k':'AM1 CHARGE=%d' % mol.charge})
mop_am1=mol.write(format='mopin',opt={'k':'AM1 CHARGE=%d' % mol.charge},overwrite=True)
print mop_am1
AM1 CHARGE=0


C    0.000000  1    0.000000  1    0.000000  1     0   0   0
C    1.347809  1    0.000000  1    0.000000  1     1   0   0
C    1.455958  1  108.591573  1    0.000000  1     2   1   0
C    1.341187  1  123.886636  1  173.842312  1     3   2   1
C    1.417910  1  123.611559  1  179.148968  1     4   3   2
C    1.380133  1  132.749293  1  147.501216  1     5   4   3
C    1.418046  1  107.757118  1  177.396528  1     6   5   4
C    1.380128  1  107.756034  1    0.000974  1     7   6   5
C    1.417955  1  132.752157  1  182.587758  1     8   7   6
C    1.340873  1  123.595469  1  212.506580  1     9   8   7
C    1.453128  1  123.873864  1  180.838426  1    10   9   8
C    1.345758  1  108.705107  1  186.162795  1    11  10   9
C    1.456045  1  108.568005  1    0.001129  1    12  11  10
C    1.341160  1  123.925407  1  173.837038  1    13  12  11
C    1.417875  1  123.590801  1  179.160826  1    14  13  12
C    1.380111  1  132.750856  1  147.504314  1    15  14  13
C    1.418055  1  107.757264  1  177.411411  1    16  15  14
C    1.380118  1  107.758682  1  359.998802  1    17  16  15
C    1.417867  1  132.747488  1  182.605651  1    18  17  16
C    1.341180  1  123.608445  1  212.488655  1    19  18  17
N    1.384486  1  130.128982  1    3.984615  1    20  19  18
N    1.375984  1  106.984219  1  358.045986  1    18  17  16
H    1.013276  1  124.796529  1  183.601560  1    22  18  17
N    1.384492  1  105.886658  1  356.333547  1    13  12  11
N    1.375931  1  120.118437  1  322.460493  1     5   4   3
H    1.013231  1  124.796047  1    0.269764  1    25   5   4
H    1.078216  1  126.223288  1  181.382343  1     1   2   3
H    1.078168  1  126.221159  1  178.616899  1     2   1   3
H    1.087321  1  119.408744  1  358.155947  1     4   3   2
H    1.082386  1  125.770018  1  357.409672  1     6   5   4
H    1.082385  1  126.473044  1  180.010348  1     7   6   5
H    1.087326  1  116.967006  1   31.530355  1     9   8   7
H    1.078247  1  125.142827  1    4.804377  1    11  10   9
H    1.078210  1  126.267419  1  178.621143  1    12  11  10
H    1.087257  1  119.417565  1  358.156375  1    14  13  12
H    1.082381  1  125.768669  1  357.422055  1    16  15  14
H    1.082383  1  126.473379  1  180.013867  1    17  16  15
H    1.087270  1  116.972738  1   31.514108  1    19  18  17
H    1.015088  1  116.042361  1   36.318173  1    21  20  19
H    1.015114  1  116.027055  1  140.870334  1    24  13  12

In [15]:
%%bash
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
/home/preps/golovin/progs/mopac/MOPAC2016.exe my.mop

          MOPAC Job: "my.mop" ended normally on May 17, 2017, at 20:45.

In [16]:
opt=pybel.readfile('mopout','my.out')
for i in opt:
    print i
    i.write(format='pdb',filename='my.pdb')
c1cc2/C=c\3/cc/c(=C\c4cc/c(=C/c5ccc(/C=c/1\[nH]2)[nH]5)/[nH]4)/[nH]3	my.out

In [20]:
import subprocess
def mopac_pr5(smiles, method, prefix):
#стоило как-то иначе назвать функцию, но не хочется все переписывать    
    mol=pybel.readstring('smi',smiles)
    mol.addh()
    mol.make3D(steps=10000)
    mop=mol.write(format='mopin',opt={'k': method+' CHARGE=%d' % mol.charge})
    mop=mol.write(format='mopin',filename = prefix+'.mop',opt={'k': method+' CHARGE=%d' % mol.charge},overwrite=True)
    mol.write(format='pdb',filename=prefix+'_before.pdb',overwrite=True)
    
    cmd1 = " export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/' "
    cmd2 = '/home/preps/golovin/progs/mopac/MOPAC2016.exe '+prefix+'.mop'
    subprocess.call(cmd1, shell=True)
    subprocess.call(cmd2, shell=True)
    opt=pybel.readfile('mopout',prefix+'.out')
    for i in opt:
        print i
        i.write(format='pdb',filename=prefix+'.pdb',overwrite=True)
In [21]:
mopac_pr5('CC(=O)C','PM6', 'acetone')
mopac_pr5('C1=CC2=CC3=CC=C(N3)C=C4C=CC(=N4)C=C5C=CC(=N5)C=C1N2','PM6', 'por_pm6')
mopac_pr5('C1=CC2=CC3=CC=C(N3)C=C4C=CC(=N4)C=C5C=CC(=N5)C=C1N2','AM1', 'por_am1')
CC(=O)C	acetone.out

c1c/c/2=C/c3cc/c(=C/c4ccc([n]4)/C=C\4/C=CC(=N4)/C=c/1\[nH]2)/[nH]3	por_pm6.out

c1c/c/2=C/c3cc/c(=C/c4ccc([n]4)/C=C\4/C=CC(=N4)/C=c/1\[nH]2)/[nH]3	por_am1.out

Получившиеся структуры:

Рис. 1. Структура порфирина, полученная MOPAC способом АМ1.

Рис. 1. Структура порфирина, полученная MOPAC способом PM6.

Из рисунок можно заметить, что правильная структура получилась в способе РМ6.

Пункт 3. Возбужденные состояния порфирина и спектр поглощения молекулы

In [22]:
%%bash

cp my.mop my_spectr.mop
echo "" >> my_spectr.mop
echo "cis c.i.=4 meci oldgeo" >> my_spectr.mop
echo "some description" >> my_spectr.mop
In [1]:
%%bash
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
/home/preps/golovin/progs/mopac/MOPAC2016.exe my_spectr.mop

          MOPAC Job: "my_spectr.mop" ended normally on May 17, 2017, at 22:13.

In [4]:
opt=pybel.readfile('mopout','my_spectr.out')
for i in opt:
    print i
    i.write(format='pdb',filename='my_spectr.pdb',overwrite=True)
c1cc2/C=c\3/cc/c(=C\c4cc/c(=C/c5ccc(/C=c/1\[nH]2)[nH]5)/[nH]4)/[nH]3	my_spectr.out

На основании значений энергий для электронных переходов, приведенных в конце файла, рассчитали длину волн, при которых происходят эти переходы ( по формуле Е=h*v)

In [6]:
f = open('my_spectr.out', 'r')
lines = f.readlines()[-21:-13]
f.close()
energies = [float(line.split()[1]) for line in lines]
print energies, 'eV'
[2.462477, 2.617337, 2.992505, 3.606955, 3.621798, 4.043092, 4.173327, 4.525836] eV

Получили из файла значения энергий, теперь рассчитаем соотвутствующие им длины волн и частоты

In [7]:
wavelengths = [1239.84193/e for e in energies]
frequencies = [e/4.135667516 for e in energies]
print 'frequncy,THz\t wavelength,nm\tenergy,eV'
for i in range(len(energies)):
        print '%.4f \t %.4f \t %.4f' % (frequencies[i]*1000, wavelengths[i], energies[i])
frequncy,THz	 wavelength,nm	energy,eV
595.4243 	 503.4938 	 2.4625
632.8693 	 473.7036 	 2.6173
723.5845 	 414.3157 	 2.9925
872.1579 	 343.7365 	 3.6070
875.7469 	 342.3277 	 3.6218
977.6153 	 306.6569 	 4.0431
1009.1060 	 297.0872 	 4.1733
1094.3423 	 273.9476 	 4.5258
In [ ]: