In [1]:
import numpy
import scipy.special
import scipy.misc
import math
import npy2cube

ДЗ выполнялось без привлечения сервера и через GUI pymol.

Общее решение уравнения Шредингера для атома водорода (волновая функция)

$\begin{equation} \psi_{nlm}(r,\theta, \varphi) = \sqrt{\frac {\left(n-l-1\right)!}{2 n \cdot\left[\left(n+l\right)!\right]^3}}\cdot\left(\frac{2}{na_0}\right)^\frac{3}{2}\cdot exp\left(-\frac{r}{na_0}\right)\cdot\left(\frac{2r}{na_0}\right)^l L^{2l+1}_{n-l-1}\left(\frac{2r}{na_0}\right)\cdot Y_{l,m}(\varphi, \theta) \end{equation}$

n,l,m - квантовые числа. n: главное квантовое число. l: орбитальное. m: магнитное.

In [176]:
def w(n,l,m,d):
    
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    # поcкольку кулоновское поле ядра атома водорода является сферически симметричным,
    #решение уравнения Шредингера ищут в сферической системе координат (r, θ, φ), и наша волновая функция (см выше) записана в сферических координатах. 
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)# переходим к сферическим координатам.
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
    phi = lambda x,y,z: numpy.arctan(y/x)
    a0 = 1.0
    
    #этот коэффициент использовать не будем, потому что с ним не получаются орбитали красивые (скорее всего придется подбирать масштаб для каждой орбитали)
    coef = lambda n,l: (numpy.sqrt((math.factorial(n-l-1))/(2.*n*(math.factorial(n+l))**3))) * ((2./n/a0)**(3./2.))#первые 2 множителя, не зависящие от радиуса и углов
    
    R = lambda r,n,l: (2.*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2.*r/n/a0)# радиальная часть волновой функции
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)# радиальная часть умножается на угловую (сферические гармоники)
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2# плотность вероятности нахождения частицы в данной точке пространства
      
    return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Боровский радиус атома водорода - 0.529 A, но это радиус ближайшей к ядру орбитали. А нас будут интересовать и следующие. При этом имеет смысл брать значение d около 2,но при таком значении не все орбитали отрисовываются. В итоге пришлось взять значение 25, хотя за ним нет никакого физического смысла.

In [177]:
d = 25.
step = float(2.*d/29)

Создадим cube-файлы для орбиталей из первых трех уровней. Вообще-то l изменяется в пределах [0,n-1], a m - [-l;l], так что мы возьмем не все орбитали.

In [178]:
for n in range(1,4):
    for l in range(0,n):
        for m in range(0,l+1):
            grid = w(n, l, m, d)
            name='%s-%s-%s' % (n,l,m)
            npy2cube.npy2cube(grid,(-d, -d, -d),(step, step, step),name+'.cube')
In [275]:
f = open ('volume_color.pml', 'w')

f.write('cmd.volume_ramp_new(\'ramp007\', [\
      -0.035, 0.00, 0.00, 1.00, 0.00, \
      -0.03,  0.00, 1.00, 1.00, 0.20, \
      -0.025, 0.00, 1.00, 0.00, 0.00, \
      0.025, 0.00, 1.00, 0.00, 0.00, \
      0.03,  0.00, 1.00, 1.00, 0.20, \
      0.035, 0.00, 0.00, 1.00, 0.00, \
    ])\n')



for n in range(1,4):
    for l in range(0,n):
        for m in range(0,l+1):
            name='%s-%s-%s' % (n,l,m)
            fname = name+'.cube'
            vol_name = name+'_v'
            f.write('load /Users/anna/Downloads/Struct_bio/pymol/%s, %s\n' % (fname, name)) 
            f.write('volume %s, %s\n' % (vol_name, name))
            f.write('volume_color %s, ramp007\n' % (vol_name)) 
f.close()

Сохраним полученные изображения и посмотрим, что получилось.

In [312]:
from IPython.display import display,Image
D=[]# list of images
for n in range(1,4):
    for l in range(n):
        for m in range(l+1):
            name='%s-%s-%s' % (n,l,m)
            fname = '/Users/anna/Downloads/Struct_bio/orbitals/'+name+'.png'
            D.append(fname)

for imageName in D:
    display(Image(imageName))

Нет уверенности в том, что такие картинки правильно показывают соотношение размеров орбиталей (s-орбитали увеличиваются при возрастании n, но, наверное, не на столько), но зато их форма вполне похожа на форму орбиталей из школьных учебников и Википедии! Для получения такой формы нужно было поменять параметры отображения volume, заданные по умолчанию в pymol (для p- и d-орбиталей рисовал только половину)T.е. создать правильный volume_ramp.

Создадим файл для orca. Для водорода попытка менять параметры заряда и мультиплетности не привела к желаемой отрисовке бОльшего количества орбиталей (только 5). Orca начинает отсчет орбиталей с нуля.

In [299]:
%%writefile h.in
! UHF SVP XYZFile
%plots Format Cube
 MO("H-0.cube",0,0);
 MO("H-1.cube",1,0);
 MO("H-2.cube",2,0);
 MO("H-3.cube",3,0);
 MO("H-4.cube",4,0);
end
* xyz 0 4
 H 0 0 0
*

Writing h.in
In [300]:
%%bash
orca '/Users/anna/Downloads/Struct_bio/pymol/h.in'

                                 *****************
                                 * 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


Your calculation utilizes the basis: Ahlrichs-VDZ
Cite in your paper:
H - Kr: A. Schaefer, H. Horn and R. Ahlrichs, J. Chem. Phys. 97, 2571 (1992).

Your calculation utilizes the basis: Ahlrichs SVPalls1+f
Cite in your paper:
Rb - Xe: A. Schaefer, C. Huber and R. Ahlrichs, J. Chem. Phys. 100, 5829 (1994).

Your calculation utilizes pol. fcns from basis: Ahlrichs polarization
Cite in your paper:
H - Kr: R. Ahlrichs and coworkers, unpublished

================================================================================
                                        WARNINGS
                       Please study these warnings very carefully!
================================================================================
Now building the actual basis set


INFO   : the flag for use of LIBINT has been found!

================================================================================
                                       INPUT FILE
================================================================================
NAME = /Users/anna/Downloads/Struct_bio/pymol/h.in
|  1> ! UHF SVP XYZFile
|  2> %plots Format Cube
|  3>  MO("H-0.cube",0,0);
|  4>  MO("H-1.cube",1,0);
|  5>  MO("H-2.cube",2,0);
|  6>  MO("H-3.cube",3,0);
|  7>  MO("H-4.cube",4,0);
|  8> end
|  9> * xyz 0 4
| 10>  H 0 0 0
| 11> **                         ****END OF INPUT****
================================================================================

                       ****************************
                       * Single Point Calculation *
                       ****************************

---------------------------------
CARTESIAN COORDINATES (ANGSTROEM)
---------------------------------
  H      0.000000    0.000000    0.000000

----------------------------
CARTESIAN COORDINATES (A.U.)
----------------------------
  NO LB      ZA    FRAG    MASS        X           Y           Z
   0 H     1.0000    0     1.008          0.000000000000000          0.000000000000000          0.000000000000000

--------------------------------
INTERNAL COORDINATES (ANGSTROEM)
--------------------------------
 H      0   0   0   0.000000     0.000     0.000

---------------------------
INTERNAL COORDINATES (A.U.)
---------------------------
 H      0   0   0   0.000000     0.000     0.000

---------------------
BASIS SET INFORMATION
---------------------
There are 1 groups of distinct atoms

 Group   1 Type H   : 4s1p contracted to 2s1p pattern {31/1}

Atom   0H    basis set group =>   1

Checking for AutoStart:
The File: /Users/anna/Downloads/Struct_bio/pymol/h.gbw exists
Trying to determine its content:
     ... Fine, the file contains calculation information
     ... Fine, the calculation information was read
     ... Fine, the file contains a basis set
     ... Fine, the basis set was read
     ... Fine, the file contains a geometry
     ... Fine, the geometry was read
     ... Fine, the file contains a set of orbitals
     ... Fine, the orbitals can be read
     => possible old guess file was deleted
     => GBW file was renamed to GES file
     => GES file is set as startup file
     => Guess is set to MORead
     ... now leaving AutoStart

------------------------------------------------------------------------------
                           ORCA GTO INTEGRAL CALCULATION
------------------------------------------------------------------------------

                         BASIS SET STATISTICS AND STARTUP INFO

 # of primitive gaussian shells          ...    5
 # of primitive gaussian functions       ...    7
 # of contracted shell                   ...    3
 # of contracted basis functions         ...    5
 Highest angular momentum                ...    1
 Maximum contraction depth               ...    3
 Integral package used                   ... LIBINT
 Integral threshhold            Thresh   ...  1.000e-10
 Primitive cut-off              TCut     ...  1.000e-11


                              INTEGRAL EVALUATION

 One electron integrals                  ... done
 Pre-screening matrix                    ... done
 Shell pair data                         ... done (   0.002 sec)

-------------------------------------------------------------------------------
                                 ORCA SCF
-------------------------------------------------------------------------------

------------
SCF SETTINGS
------------
Hamiltonian:
 Ab initio Hamiltonian  Method          .... Hartree-Fock(GTOs)


General Settings:
 Integral files         IntName         .... /Users/anna/Downloads/Struct_bio/pymol/h
 Hartree-Fock type      HFTyp           .... UHF
 Total Charge           Charge          ....    0
 Multiplicity           Mult            ....    4
 Number of Electrons    NEL             ....    1
 Basis Dimension        Dim             ....    5
 Nuclear Repulsion      ENuc            ....      0.0000000000 Eh

Convergence Acceleration:
 DIIS                   CNVDIIS         .... on
   Start iteration      DIISMaxIt       ....    12
   Startup error        DIISStart       ....  0.200000
   # of expansion vecs  DIISMaxEq       ....     5
   Bias factor          DIISBfac        ....   1.050
   Max. coefficient     DIISMaxC        ....  10.000
 Newton-Raphson         CNVNR           .... off
 SOSCF                  CNVSOSCF        .... off
 Level Shifting         CNVShift        .... on
   Level shift para.    LevelShift      ....    0.2500
   Turn off err/grad.   ShiftErr        ....    0.0010
 Zerner damping         CNVZerner       .... off
 Static damping         CNVDamp         .... on
   Fraction old density DampFac         ....    0.7000
   Max. Damping (<1)    DampMax         ....    0.9800
   Min. Damping (>=0)   DampMin         ....    0.0000
   Turn off err/grad.   DampErr         ....    0.1000
 Fernandez-Rico         CNVRico         .... off

SCF Procedure:
 Maximum # iterations   MaxIter         ....   125
 SCF integral mode      SCFMode         .... Direct
   Integral package                     .... LIBINT
 Reset frequeny         DirectResetFreq ....    20
 Integral Threshold     Thresh          ....  1.000e-10 Eh
 Primitive CutOff       TCut            ....  1.000e-11 Eh

Convergence Tolerance:
 Convergence Check Mode ConvCheckMode   .... Total+1el-Energy
 Energy Change          TolE            ....  1.000e-06 Eh
 1-El. energy change                    ....  1.000e-03 Eh
 DIIS Error             TolErr          ....  1.000e-06


Diagonalization of the overlap matrix:
Smallest eigenvalue                        ... 3.152e-01
Time for diagonalization                   ...    0.000 sec
Threshold for overlap eigenvalues          ... 1.000e-08
Number of eigenvalues below threshold      ... 0
Time for construction of square roots      ...    0.001 sec
Total time needed                          ...    0.001 sec

---------------------
INITIAL GUESS: MOREAD
---------------------
Guess MOs are being read from file: /Users/anna/Downloads/Struct_bio/pymol/h.ges
Input Geometry matches current geometry (good)
Atom 0: N(Shells)= 3 and 8 - projection required
Input basis set needs projection
NBasis(read)= 18
NBasis(calc)= 5
Basis set overlap was calculated
Using GuessMode=FMatrix for the initial guess
C*T matrix was formed
Effective Fockian was calculated
Effective Fockian was diagonalized
Initial density was built
C*T matrix was formed
Effective Fockian was calculated
Effective Fockian was diagonalized
Initial density was built
                      ------------------
                      INITIAL GUESS DONE (   0.0 sec)
                      ------------------
--------------
SCF ITERATIONS
--------------
ITER       Energy         Delta-E        Max-DP      RMS-DP      [F,P]     Damp
               ***  Starting incremental Fock matrix formation  ***
  0      1.5700552820   0.000000000000 0.46782962  0.07191201  0.2296337 0.7000
  1      1.2011895997  -0.368865682323 0.32748074  0.05033841  0.1646005 0.7000
                               ***Turning on DIIS***
  2      0.9020346253  -0.299154974418 0.22923651  0.03523689  0.1171103 0.7000
  3      0.6242574267  -0.277777198548 0.53488520  0.08221940  0.0829032 0.0000
                            ***DIIS convergence achieved***

               *****************************************************
               *                     SUCCESS                       *
               *           SCF CONVERGED AFTER   4 CYCLES          *
               *****************************************************


----------------
TOTAL SCF ENERGY
----------------

Total Energy       :            0.07286251 Eh               1.98269 eV

Components:
Nuclear Repulsion  :            0.00000000 Eh               0.00000 eV
Electronic Energy  :            0.07286251 Eh               1.98269 eV

One Electron Energy:           -0.31757803 Eh              -8.64174 eV
Two Electron Energy:            0.39044053 Eh              10.62443 eV

Virial components:
Potential Energy   :           -1.57795800 Eh             -42.93842 eV
Kinetic Energy     :            1.65082051 Eh              44.92111 eV
Virial Ratio       :            0.95586285


---------------
SCF CONVERGENCE
---------------

  Last Energy change         ...   -5.5139e-01  Tolerance :   1.0000e-06
  Last MAX-Density change    ...    7.6067e-09  Tolerance :   1.0000e-05
  Last RMS-Density change    ...    1.7477e-09  Tolerance :   1.0000e-06
  Last DIIS Error            ...    8.2837e-09  Tolerance :   1.0000e-06

             **** THE GBW FILE WAS UPDATED (/Users/anna/Downloads/Struct_bio/pymol/h.gbw) ****
             **** DENSITY FILE WAS UPDATED (/Users/anna/Downloads/Struct_bio/pymol/h.scfp.tmp) ****
             **** ENERGY FILE WAS UPDATED (/Users/anna/Downloads/Struct_bio/pymol/h.en.tmp) ****
----------------------
UHF SPIN CONTAMINATION
----------------------

Expectation value of <S**2>     :     2.000000
Ideal value S*(S+1) for S=1.0   :     2.000000
Deviation                       :     0.000000

----------------
ORBITAL ENERGIES
----------------
                 SPIN UP ORBITALS
  NO   OCC          E(Eh)            E(eV) 
   0   1.0000      -0.108838        -2.9616 
   1   1.0000       0.572141        15.5687 
   2   0.0000       2.100372        57.1540 
   3   0.0000       2.100372        57.1540 
   4   0.0000       2.100372        57.1540 

                 SPIN DOWN ORBITALS
  NO   OCC          E(Eh)            E(eV) 
   0   0.0000       0.487405        13.2630 
   1   0.0000       1.318415        35.8759 
   2   0.0000       2.270122        61.7732 
   3   0.0000       2.270122        61.7732 
   4   0.0000       2.270122        61.7732 

                    ********************************
                    * MULLIKEN POPULATION ANALYSIS *
                    ********************************

      **** WARNING: MULLIKEN FINDS    2.0000000 ELECTRONS INSTEAD OF 1 ****
--------------------------------------------
MULLIKEN ATOMIC CHARGES AND SPIN POPULATIONS
--------------------------------------------
   0 H :   -1.000000    2.000000
Sum of atomic charges         :   -1.0000000
Sum of atomic spin populations:    2.0000000

-----------------------------------------------------
MULLIKEN REDUCED ORBITAL CHARGES AND SPIN POPULATIONS
-----------------------------------------------------
CHARGE
  0 H s       :     2.000000  s :     2.000000
      pz      :     0.000000  p :     0.000000
      px      :     0.000000
      py      :     0.000000

SPIN
  0 H s       :     2.000000  s :     2.000000
      pz      :     0.000000  p :     0.000000
      px      :     0.000000
      py      :     0.000000


                     *******************************
                     * LOEWDIN POPULATION ANALYSIS *
                     *******************************

      **** WARNING: LOEWDIN FINDS    2.0000000 ELECTRONS INSTEAD OF 1 ****
-------------------------------------------
LOEWDIN ATOMIC CHARGES AND SPIN POPULATIONS
-------------------------------------------
   0 H :   -1.000000    2.000000

----------------------------------------------------
LOEWDIN REDUCED ORBITAL CHARGES AND SPIN POPULATIONS
----------------------------------------------------
CHARGE
  0 H s       :     2.000000  s :     2.000000
      pz      :     0.000000  p :     0.000000
      px      :     0.000000
      py      :     0.000000

SPIN
  0 H s       :     2.000000  s :     2.000000
      pz      :     0.000000  p :     0.000000
      px      :     0.000000
      py      :     0.000000


                      *****************************
                      * MAYER POPULATION ANALYSIS *
                      *****************************

  NA   - Mulliken gross atomic population
  ZA   - Total nuclear charge
  QA   - Mulliken gross atomic charge
  VA   - Mayer's total valence
  BVA  - Mayer's bonded valence
  FA   - Mayer's free valence

  ATOM       NA         ZA         QA         VA         BVA        FA
  0 H      2.0000     1.0000    -1.0000     2.0000     0.0000     2.0000

  Mayer bond orders larger than 0.1


-----------------------------------------------
ATOM BASIS FOR ELEMENT H 
-----------------------------------------------
  NSH[1]   = 3;
  res=GAUSS_InitGTOSTO(BG,BS,1,NSH[1]);
  // Basis function for L=s
  (*BG)[ 1][ 0].l    = ((*BS)[ 1][ 0].l=0);
  (*BG)[ 1][ 0].ng   = 4;
  (*BG)[ 1][ 0].a[ 0]=     13.01070100; (*BG)[ 1][ 0].d[ 0]= 0.096096678877;
  (*BG)[ 1][ 0].a[ 1]=      1.96225720; (*BG)[ 1][ 0].d[ 1]= 0.163022191701;
  (*BG)[ 1][ 0].a[ 2]=      0.44453796; (*BG)[ 1][ 0].d[ 2]= 0.185592186247;
  (*BG)[ 1][ 0].a[ 3]=      0.12194962; (*BG)[ 1][ 0].d[ 3]= 0.073701452542;
  // Basis function for L=s
  (*BG)[ 1][ 1].l    = ((*BS)[ 1][ 1].l=0);
  (*BG)[ 1][ 1].ng   = 4;
  (*BG)[ 1][ 1].a[ 0]=     13.01070100; (*BG)[ 1][ 1].d[ 0]= 0.202726593388;
  (*BG)[ 1][ 1].a[ 1]=      1.96225720; (*BG)[ 1][ 1].d[ 1]= 0.343913379281;
  (*BG)[ 1][ 1].a[ 2]=      0.44453796; (*BG)[ 1][ 1].d[ 2]= 0.391527283951;
  (*BG)[ 1][ 1].a[ 3]=      0.12194962; (*BG)[ 1][ 1].d[ 3]=-0.187888280974;
  // Basis function for L=p
  (*BG)[ 1][ 2].l    = ((*BS)[ 1][ 2].l=1);
  (*BG)[ 1][ 2].ng   = 3;
  (*BG)[ 1][ 2].a[ 0]=      0.80000000; (*BG)[ 1][ 2].d[ 0]= 0.856608271641;
  (*BG)[ 1][ 2].a[ 1]=      0.80000000; (*BG)[ 1][ 2].d[ 1]= 0.476128394479;
  (*BG)[ 1][ 2].a[ 2]=      0.80000000; (*BG)[ 1][ 2].d[ 2]= 0.450102341869;
 newgto H 
 S 4
    1         13.010701000000         0.019682160277
    2          1.962257200000         0.137965241943
    3          0.444537960000         0.478319356737
    4          0.121949620000         0.501107169599
 S 4
    1         13.010701000000         0.041521698254
    2          1.962257200000         0.291052966991
    3          0.444537960000         1.009067689709
    4          0.121949620000        -1.277480448918
 P 3
    1          0.800000000000         0.794291092301
    2          0.800000000000         0.441490649864
    3          0.800000000000         0.417357959999
 end
-------------------------------------------
RADIAL EXPECTATION VALUES <R**-3> TO <R**3>
-------------------------------------------
   0 :     0.000000     1.965412     0.998557     1.497110     2.972853     7.305452
   1 :     0.000000     2.769720     0.969842     2.227011     6.779696    23.235533
   2 :     1.522453     1.066667     0.951533     1.189416     1.562500     2.230155
   3 :     1.522453     1.066667     0.951533     1.189416     1.562500     2.230155
   4 :     1.522453     1.066667     0.951533     1.189416     1.562500     2.230155
-------
TIMINGS
-------

Total SCF time: 0 days 0 hours 0 min 1 sec 

Total time                  ....       1.112 sec
Sum of individual times     ....       1.108 sec  ( 99.7%)

Fock matrix formation       ....       1.093 sec  ( 98.3%)
Diagonalization             ....       0.001 sec  (  0.0%)
Density matrix formation    ....       0.000 sec  (  0.0%)
Population analysis         ....       0.008 sec  (  0.7%)
Initial guess               ....       0.003 sec  (  0.3%)
Orbital Transformation      ....       0.000 sec  (  0.0%)
Orbital Orthonormalization  ....       0.000 sec  (  0.0%)
DIIS solution               ....       0.004 sec  (  0.4%)

-------------------------   --------------------
FINAL SINGLE POINT ENERGY         0.072862505694
-------------------------   --------------------


---------------
PLOT GENERATION
---------------
choosing x-range =    -7.000000 ..     7.000000
choosing y-range =    -7.000000 ..     7.000000
choosing z-range =    -7.000000 ..     7.000000

GBW-File       ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
PlotType       ... MO-PLOT
MO/Operator    ... 0 0
Output file    ... H-0.cube
Format         ... Gaussian-Cube
Resolution     ... 40 40 40
Boundaries     ...    -7.000000     7.000000 (x direction)
                      -7.000000     7.000000 (y direction)
                      -7.000000     7.000000 (z direction)
choosing x-range =    -7.000000 ..     7.000000
choosing y-range =    -7.000000 ..     7.000000
choosing z-range =    -7.000000 ..     7.000000

GBW-File       ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
PlotType       ... MO-PLOT
MO/Operator    ... 1 0
Output file    ... H-1.cube
Format         ... Gaussian-Cube
Resolution     ... 40 40 40
Boundaries     ...    -7.000000     7.000000 (x direction)
                      -7.000000     7.000000 (y direction)
                      -7.000000     7.000000 (z direction)
choosing x-range =    -7.000000 ..     7.000000
choosing y-range =    -7.000000 ..     7.000000
choosing z-range =    -7.000000 ..     7.000000

GBW-File       ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
PlotType       ... MO-PLOT
MO/Operator    ... 2 0
Output file    ... H-2.cube
Format         ... Gaussian-Cube
Resolution     ... 40 40 40
Boundaries     ...    -7.000000     7.000000 (x direction)
                      -7.000000     7.000000 (y direction)
                      -7.000000     7.000000 (z direction)
choosing x-range =    -7.000000 ..     7.000000
choosing y-range =    -7.000000 ..     7.000000
choosing z-range =    -7.000000 ..     7.000000

GBW-File       ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
PlotType       ... MO-PLOT
MO/Operator    ... 3 0
Output file    ... H-3.cube
Format         ... Gaussian-Cube
Resolution     ... 40 40 40
Boundaries     ...    -7.000000     7.000000 (x direction)
                      -7.000000     7.000000 (y direction)
                      -7.000000     7.000000 (z direction)
choosing x-range =    -7.000000 ..     7.000000
choosing y-range =    -7.000000 ..     7.000000
choosing z-range =    -7.000000 ..     7.000000

GBW-File       ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
PlotType       ... MO-PLOT
MO/Operator    ... 4 0
Output file    ... H-4.cube
Format         ... Gaussian-Cube
Resolution     ... 40 40 40
Boundaries     ...    -7.000000     7.000000 (x direction)
                      -7.000000     7.000000 (y direction)
                      -7.000000     7.000000 (z direction)

                            ***************************************
                            *     ORCA property calculations      *
                            ***************************************

                                    ---------------------
                                    Active property flags
                                    ---------------------
   (+) Dipole Moment


------------------------------------------------------------------------------
                       ORCA ELECTRIC PROPERTIES CALCULATION
------------------------------------------------------------------------------

Dipole Moment Calculation                       ... on
Quadrupole Moment Calculation                   ... off
Polarizability Calculation                      ... off
GBWName                                         ... /Users/anna/Downloads/Struct_bio/pymol/h.gbw
Electron density file                           ... /Users/anna/Downloads/Struct_bio/pymol/h.scfp.tmp

-------------
DIPOLE MOMENT
-------------
                                X             Y             Z
Electronic contribution:     -0.00000       0.00000       0.00000
Nuclear contribution   :      0.00000       0.00000       0.00000
                        -----------------------------------------
Total Dipole Moment    :     -0.00000       0.00000       0.00000
                        -----------------------------------------
Magnitude (a.u.)       :      0.00000
Magnitude (Debye)      :      0.00000


Timings for individual modules:

Sum of individual times         ...        3.650 sec (=   0.061 min)
GTO integral calculation        ...        1.081 sec (=   0.018 min)  29.6 %
SCF iterations                  ...        1.620 sec (=   0.027 min)  44.4 %
Orbital/Density plot generation ...        0.949 sec (=   0.016 min)  26.0 %
                             ****ORCA TERMINATED NORMALLY****
TOTAL RUN TIME: 0 days 0 hours 0 minutes 4 seconds 528 msec
In [301]:
f2 = open ('volume_color_orca.pml', 'w')

f2.write('cmd.volume_ramp_new(\'ramp008\', [\
      -0.035, 0.00, 0.00, 1.00, 0.00, \
      -0.03,  0.00, 1.00, 1.00, 0.20, \
      -0.025, 0.00, 1.00, 0.00, 0.00, \
      0.025, 0.00, 1.00, 0.00, 0.00, \
      0.03,  0.00, 1.00, 1.00, 0.20, \
      0.035, 0.00, 0.00, 1.00, 0.00, \
    ])\n')


for i in range(0,5):
    
            name="H-"+'%s' % (i)
            fname = name+'.cube'
            vol_name = name+'_V'
            f2.write('load /Users/anna/Downloads/Struct_bio/pymol/%s, %s\n' % (fname, name)) 
            f2.write('volume %s, %s\n' % (vol_name, name))
            f2.write('volume_color %s, ramp008\n' % (vol_name)) 
f2.close()
In [307]:
display(Image('/Users/anna/Downloads/Struct_bio/orbitals/H-0.png'))
display(Image('/Users/anna/Downloads/Struct_bio/orbitals/H-1.png'))
display(Image('/Users/anna/Downloads/Struct_bio/orbitals/H-2.png'))
display(Image('/Users/anna/Downloads/Struct_bio/orbitals/H-3.png'))

Полученные орбитали похожи на те, что мы рассчитывали сами!

Также были построены орбитали для атома с бОльшим количеством электронов - для Al, чтобы посмотреть, что происходит на более высоком уровне. Не все орбитали на уровне выше 2, полученные для Al, похожи на рассчитанные орбитали H.