1) модуль pybel

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

Неправильная структура: не плоская и водороды в центре.

2). MOPAC способ PM6/АМ1

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

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

In [6]:
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 [8]:
%%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 25, 2017, at 18:13.

In [9]:
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 [2]:
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 [11]:
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

Методом AM1 не получился плоский порфирин, а PM6 получился. На рисунке - наложение структур, беленький - не плоский.

In [1]:
%%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 [2]:
%%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 25, 2017, at 20:32.

In [5]:
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

In [24]:
from IPython.core.display import HTML 
In [8]:
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 [9]:
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 [11]:
import pybel

mol=pybel.readfile("pdb", "td.pdb").next()
mop=mol.write(format='mop',filename='td.mop',opt={'k':'PM6 CHARGE=%d' % mol.charge},overwrite=True)
In [12]:
%%bash
export MOPAC_LICENSE='/home/preps/golovin/progs/mopac/'
/home/preps/golovin/progs/mopac/MOPAC2016.exe td.mop

          MOPAC Job: "td.mop" ended normally on May 25, 2017, at 20:43.

In [13]:
opt=pybel.readfile('mopout','td.out')
for i in opt:
    print i
    i.write(format='pdb',filename='td_opt.pdb',overwrite=True)
N1[C@H]2[C@@](C)(C(=O)NC1=O)[C@@]1([C@H]2NC(=O)NC1=O)C	td.out

Исходная и оптимизированная (розовенькая) структуры.

In [14]:
%%bash
cp td.mop td2.mop
#change line PM6 CHARGE=0 to PM6 CHARGE=+2
In [15]:
! /home/preps/golovin/progs/mopac/MOPAC2016.exe td2.mop

          MOPAC Job: "td2.mop" ended normally on May 25, 2017, at 20:44.

In [16]:
opt=pybel.readfile('mopout','td2.out')
for i in opt:
    print i
    i.write(format='pdb',filename='td2_opt.pdb',overwrite=True)
N1[C@H]2[C@@](C)(C(=O)NC1=O)[C@@]1([C@H]2NC(=O)NC1=O)C	td2.out

Два тимина, соединены все равно двумя связями.

In [17]:
mol=pybel.readfile("pdb", "td2_opt.pdb").next()
mop=mol.write(format='mop',filename='td2_opt.mop',opt={'k':'PM6 CHARGE=%d' % mol.charge},overwrite=True)
In [18]:
! /home/preps/golovin/progs/mopac/MOPAC2016.exe td2_opt.mop

          MOPAC Job: "td2_opt.mop" ended normally on May 25, 2017, at 20:45.

In [19]:
opt=pybel.readfile('mopout','td2_opt.out')
for i in opt:
    print i
    i.write(format='pdb',filename='td2_last.pdb',overwrite=True)
N1[C@H]2[C@@](C)(C(=O)NC1=O)[C@@]1([C@H]2NC(=O)NC1=O)C	td2_opt.out

Оптимизировали их:

In [20]:
%%bash
echo dimer
grep 'TOTAL ENERGY' td.out
echo semi-dimer
grep 'TOTAL ENERGY' td2.out
echo two thymins
grep 'TOTAL ENERGY' td2_opt.out
dimer
          TOTAL ENERGY            =      -3273.57488 EV
semi-dimer
          TOTAL ENERGY            =      -3273.57488 EV
two thymins
          TOTAL ENERGY            =      -3273.57524 EV

Ab initio вычисления: Gamess US

In [5]:
mopac_pr5('C1=CC=C2C=CC=CC2=C1','PM6','na')
mopac_pr5('C1=CC=C2C=CC=C2C=C1','PM6','az')
c1ccc2ccccc2c1	na.out

c1ccc2ccccc2c1	az.out

Азален до и после оптимизации (где какой и так понятно):

Нафтален до (светло-розовый) и после оптимизации:

In [6]:
mopac_pr5('c1ccc2ccccc2c1','PM6','na_ar')
mopac_pr5('c1ccc2cccc2cc1','PM6','az_ar')
c1ccc2ccccc2c1	na_ar.out

c1ccc2ccccc2c1	az_ar.out

Ароматический азулен до и после оптимизации: Ароматический нафтален до и после оптимизации (фиолетовый):

In [7]:
opt=pybel.readfile('mopout','na_ar.out')
for i in opt:
    print i
    i.write(format='orcainp',filename='orca.inp',overwrite=True)
c1ccc2ccccc2c1	na_ar.out

In [12]:
%%bash
/srv/databases/orca/orca orca.inp
                                 *****************
                                 * O   R   C   A *
                                 *****************

           --- An Ab Initio, DFT and Semiempirical electronic structure package ---

                  #######################################################
                  #                        -***-                        #
                  #  Department of molecular theory and spectroscopy    #
                  #              Directorship: Frank Neese              #
                  # Max Planck Institute for Chemical Energy Conversion #
                  #                  D-45470 Muelheim/Ruhr              #
                  #                       Germany                       #
                  #                                                     #
                  #                  All rights reserved                #
                  #                        -***-                        #
                  #######################################################


                         Program Version 3.0.3 - RELEASE   -


 With contributions from (in alphabetic order):
   Ute Becker             : Parallelization
   Dmytro Bykov           : SCF Hessian
   Dmitry Ganyushin       : Spin-Orbit,Spin-Spin,Magnetic field MRCI
   Andreas Hansen         : Spin unrestricted coupled pair/coupled cluster methods
   Dimitrios Liakos       : Extrapolation schemes; parallel MDCI
   Robert Izsak           : Overlap fitted RIJCOSX, COSX-SCS-MP3
   Christian Kollmar      : KDIIS, OOCD, Brueckner-CCSD(T), CCSD density
   Simone Kossmann        : Meta GGA functionals, TD-DFT gradient, OOMP2, MP2 Hessian
   Taras Petrenko         : DFT Hessian,TD-DFT gradient, ASA and ECA modules, normal mode analysis, Resonance Raman, ABS, FL, XAS/XES, NRVS
   Christoph Reimann      : Effective Core Potentials
   Michael Roemelt        : Restricted open shell CIS
   Christoph Riplinger    : Improved optimizer, TS searches, QM/MM, DLPNO-CCSD
   Barbara Sandhoefer     : DKH picture change effects
   Igor Schapiro          : Molecular dynamics
   Kantharuban Sivalingam : CASSCF convergence, NEVPT2
   Boris Wezisla          : Elementary symmetry handling
   Frank Wennmohs         : Technical directorship


 We gratefully acknowledge several colleagues who have allowed us to
 interface, adapt or use parts of their codes:
   Stefan Grimme, W. Hujo, H. Kruse, T. Risthaus : VdW corrections, initial TS optimization,
                                                   DFT functionals, gCP
   Ed Valeev                                     : LibInt (2-el integral package), F12 methods
   Garnet Chan, S. Sharma, R. Olivares           : DMRG
   Ulf Ekstrom                                   : XCFun DFT Library
   Mihaly Kallay                                 : mrcc  (arbitrary order and MRCC methods)
   Andreas Klamt, Michael Diedenhofen            : otool_cosmo (COSMO solvation model)
   Frank Weinhold                                : gennbo (NPA and NBO analysis)
   Christopher J. Cramer and Donald G. Truhlar   : smd solvation model


 Your calculation uses the libint2 library for the computation of 2-el integrals
 For citations please refer to: http://libint.valeyev.net

 This ORCA versions uses:
   CBLAS   interface :  Fast vector & matrix operations
   LAPACKE interface :  Fast linear algebra routines
   SCALAPACK package :  Parallel linear algebra routines



            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                                 INPUT ERROR
            UNRECOGNIZED OR DUPLICATED KEYWORD(S) IN SIMPLE INPUT LINE
              INSERT INLINE COMMANDS HERE  
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


Умываю руки.