Вычисление атомных орбиталей водорода

In [19]:
import numpy
import scipy.special
import scipy.misc
from IPython.display import display,Image

import os, sys

Полезная функция:

In [27]:
import npy2cube

Выглядит так:

def npy2cube(grid, start, step, cube_f):

    '''
    PARAMETERS:
        grid : numpy array
            3-dimentional array, containing grid data
        start : tuple
            format: (x, y, z), coordinates of cube start point
        step: tuple
            format: (x, y, z), step size on 3 axes 
        cube_f: string
             name of output .cube file

    RETURNS:
        void    
    '''

    cube_string = ""
    bohr_to_angs = 0.529177210859 #const
    start = map(lambda x: x/bohr_to_angs, start)
    step = map(lambda x: x/bohr_to_angs, step)

    with open(cube_f, "w") as cube_file:

        ###HEADER###
        cube_file.write(" CPMD CUBE FILE.\nOUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n")
        cube_file.write("    1 %f %f %f\n" %(start[0], start[1], start[2]))
        cube_file.write("  %i    %f   0.000000     0.000000\n" %(grid.shape[0], step[0]))
        cube_file.write("  %i    0.000000    %f    0.000000\n" %(grid.shape[1], step[1]))
        cube_file.write("  %i    0.000000    0.000000    %f\n" %(grid.shape[2], step[2]))
        cube_file.write("    1    0.000000    %f %f %f\n" %(start[0], start[1], start[2]))

        ###DATA###
        i = 0
        for x in range(grid.shape[0]):
            for y in range(grid.shape[1]):
                for z in range(grid.shape[2]):
                    if i < 5:
                        cube_file.write("%f " %(float(grid[x, y, z])))
                        i += 1
                    elif i == 5:
                        cube_file.write("%f\n" %(float(grid[x, y, z])))
                        i = 0
    return 0    


if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description = "Convert .npy to .cube file")
    parser.add_argument("-npy", help = "Input .npy filename")
    parser.add_argument("-cube", help = "Output .cube filename")
    parser.add_argument("-start", type = float, nargs = 3, help = "Coordinates of cube start point")
    parser.add_argument("-step", type = float, nargs = 3, help = "Step between nodes in grid on each axe, Angstremes")
    args = parser.parse_args()
    grid = np.load(args.npy)
    npy2grid(grid, args.start, args.step, args.cube)

z = Image(filename='sc.png') display(z)

In [28]:
# https://en.wikipedia.org/wiki/Wave_function

def w(n,l,m,d):
    
    # Grid
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    
    # Spherical coordinates.    
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)      # radial distance    
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))    # the polar angle
    phi = lambda x,y,z: numpy.arctan(y/x)             # the azimuth angle
    
    # Boron radius
    a0 = 1.
    
    # R - radial functions,
    # genlaguerre - L2ℓ + 1 n − ℓ − 1 are the generalized Laguerre polynomials of degree n − ℓ − 1,
    # n = 1, 2, ... - the principal quantum number,
    # ℓ = 0, 1, ... n − 1 - the azimuthal quantum number, 
    # m = −ℓ, −ℓ + 1, ..., ℓ − 1, ℓ - the magnetic quantum number
    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)

    # Coefficient
    sq = lambda n,l: numpy.sqrt((factorial(n-l-1)/(2*n*factorial(n+l)))\
                                        * (2/n/a0)**3)    
    # WF - wave function,
    # sph_harm - Ym ℓ(θ, φ) are spherical harmonics of degree ℓ and order m
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
    # Probability density
    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)
In [29]:
# Функция w выдает трехмерный массив из 30*30*30 элементов с шагом от -d до +d (или grid).

# Функции для первых трех уровней
step = 1

for n in range(1,4):
    d = 10 * n
    for l in range(0,n):
        for m in range(0,l+1,1):
            grid= w(n,l,m,d)
            name='%s-%s-%s' % (n, l, m)
            
            #  Запись в cube-файлы: cube_file.write("%f " %(float(grid[x, y, z])))
            npy2cube.npy2cube(grid, (-d,-d,-d),(step, step, step), name + '.cube')

npy2cube.py:45: ComplexWarning: Casting complex values to real discards the imaginary part
npy2cube.py:48: ComplexWarning: Casting complex values to real discards the imaginary part
In [30]:
x = Image(filename='h1.jpg')
display(x)
In [31]:
step = 1    
    
# Функции для более высоких уровней    
for n in range(15,20):
    d = 10 * n
    for l in range(0,n):
        for m in range(0,l+1,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 [32]:
y = Image(filename='h2.jpg')
display(y)
In [16]:
h = '''! UHF SVP XYZFile
%plots Format Cube
 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
*'''

with open('h.inp', 'w+') as wf:
    wf.write(h)
In [17]:
! cat h.inp

! UHF SVP XYZFile
%plots Format Cube
 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
*
In [10]:
%%bash
export PATH=${PATH}:/home/domain/data/prog/orca_4_0_1_2_linux_x86-64_openmpi202
orca h.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 4.0.1.2 - RELEASE -


 With contributions from (in alphabetic order):
   Daniel Aravena         : Magnetic Properties
   Michael Atanasov       : Ab Initio Ligand Field Theory
   Ute Becker             : Parallelization
   Martin Brehm           : Molecular dynamics
   Dmytro Bykov           : SCF Hessian
   Vijay G. Chilkuri      : MRCI spin determinant printing
   Dipayan Datta          : RHF DLPNO-CCSD density
   Achintya Kumar Dutta   : EOM-CC, STEOM-CC
   Dmitry Ganyushin       : Spin-Orbit,Spin-Spin,Magnetic field MRCI
   Yang Guo               : DLPNO-NEVPT2, CIM, IAO-localization
   Andreas Hansen         : Spin unrestricted coupled pair/coupled cluster methods
   Lee Huntington         : MR-EOM, pCC
   Robert Izsak           : Overlap fitted RIJCOSX, COSX-SCS-MP3, EOM
   Christian Kollmar      : KDIIS, OOCD, Brueckner-CCSD(T), CCSD density
   Simone Kossmann        : Meta GGA functionals, TD-DFT gradient, OOMP2, MP2 Hessian
   Martin Krupicka        : AUTO-CI
   Dagmar Lenk            : GEPOL surface
   Dimitrios Liakos       : Extrapolation schemes; parallel MDCI
   Dimitrios Manganas     : ROCIS; embedding schemes
   Dimitrios Pantazis     : SARC Basis sets
   Taras Petrenko         : DFT Hessian,TD-DFT gradient, ASA, ECA, R-Raman, ABS, FL, XAS/XES, NRVS
   Peter Pinski           : DLPNO-MP2
   Christoph Reimann      : Effective Core Potentials
   Marius Retegan         : Local ZFS, SOC
   Christoph Riplinger    : Optimizer, TS searches, QM/MM, DLPNO-CCSD(T), (RO)-DLPNO pert. Triples
   Tobias Risthaus        : Range-separated hybrids, TD-DFT gradient, RPA, STAB
   Michael Roemelt        : Restricted open shell CIS
   Masaaki Saitow         : Open-shell DLPNO
   Barbara Sandhoefer     : DKH picture change effects
   Kantharuban Sivalingam : CASSCF convergence, NEVPT2, FIC-MRCI
   Georgi Stoychev        : AutoAux
   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,             : VdW corrections, initial TS optimization,
                  C. Bannwarth                     DFT functionals, gCP, sTDA/sTD-DF
   Ed Valeev                                     : LibInt (2-el integral package), F12 methods
   Garnet Chan, S. Sharma, J. Yang, 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)
   Jiri Pittner, Ondrej Demel                    : Mk-CCSD
   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: SVP 

================================================================================
                                        WARNINGS
                       Please study these warnings very carefully!
================================================================================

Warning: TCutStore was < 0. Adjusted to Thresh (uncritical)

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

================================================================================
                                       INPUT FILE
================================================================================
NAME = h.inp
|  1> ! UHF SVP XYZFile
|  2> %plots Format Cube
|  3>  MO("H-1.cube",1,0);
|  4>  MO("H-2.cube",2,0);
|  5>  MO("H-3.cube",3,0);
|  6>  MO("H-4.cube",4,0);
|  7> end
|  8> * xyz 0 4
|  9>  H 0 0 0
| 10> **                         ****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.000000    0.000000    0.000000

--------------------------------
INTERNAL COORDINATES (ANGSTROEM)
--------------------------------
 H      0   0   0     0.000000000000     0.00000000     0.00000000

---------------------------
INTERNAL COORDINATES (A.U.)
---------------------------
 H      0   0   0     0.000000000000     0.00000000     0.00000000

---------------------
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
------------------------------------------------------------------------------
                           ORCA GTO INTEGRAL CALCULATION
------------------------------------------------------------------------------

                         BASIS SET STATISTICS AND STARTUP INFO

 # of primitive gaussian shells          ...    5
 # of primitive gaussian functions       ...    7
 # of contracted shells                  ...    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 
 Pre-screening matrix                    ... done
 Shell pair data                         ... done (   0.000 sec)

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

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


General Settings:
 Integral files         IntName         .... 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
 Convergence forced     ConvForced      .... 0
 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.017 sec
Threshold for overlap eigenvalues          ... 1.000e-08
Number of eigenvalues below threshold      ... 0
Time for construction of square roots      ...    0.014 sec
Total time needed                          ...    0.030 sec

-------------------
DFT GRID GENERATION
-------------------

General Integration Accuracy     IntAcc      ...  4.340
Radial Grid Type                 RadialGrid  ... Gauss-Chebyshev
Angular Grid (max. acc.)         AngularGrid ... Lebedev-110
Angular grid pruning method      GridPruning ... 3 (G Style)
Weight generation scheme         WeightScheme... Becke
Basis function cutoff            BFCut       ...    1.0000e-10
Integration weight cutoff        WCut        ...    1.0000e-14
Grids for H and He will be reduced by one unit

# of grid points (after initial pruning)     ...    794 (   0.0 sec)
# of grid points (after weights+screening)   ...    794 (   0.0 sec)
Grid point division into batches done        ...    0.0 sec
Reduced shell lists constructed in    0.0 sec

Total number of grid points                  ...      794
Total number of batches                      ...       13
Average number of points per batch           ...       61
Average number of grid points per atom       ...      794
Average number of shells per batch           ...     2.79 (92.86%)
Average number of basis functions per batch  ...     4.64 (92.86%)
Average number of large shells per batch     ...     2.79 (100.00%)
Average number of large basis fcns per batch ...     4.64 (100.00%)
Maximum spatial batch extension              ...  17.62, 21.59, 21.59 au
Average spatial batch extension              ...   6.10,  8.93,  9.49 au

Time for grid setup =    0.007 sec

------------------------------
INITIAL GUESS: MODEL POTENTIAL
------------------------------
Loading Hartree-Fock densities                     ... done
Calculating cut-offs                               ... done
Setting up the integral package                    ... done
Initializing the effective Hamiltonian             ... done
Starting the Coulomb interaction                   ... done (   0.0 sec)
Reading the grid                                   ... done
Mapping shells                                     ... done
Starting the XC term evaluation                    ... done (   0.0 sec)
Transforming the Hamiltonian                       ... done (   0.0 sec)
Diagonalizing the Hamiltonian                      ... done (   0.0 sec)
Back transforming the eigenvectors                 ... done (   0.0 sec)
Now organizing SCF variables                       ... done
                      ------------------
                      INITIAL GUESS DONE (   0.2 sec)
                      ------------------
--------------
SCF ITERATIONS
--------------
ITER       Energy         Delta-E        Max-DP      RMS-DP      [F,P]     Damp
               ***  Starting incremental Fock matrix formation  ***
  0      0.0728625057   0.000000000000 0.00000000  0.00000000  0.0000000 0.7000
                 **** Energy Check signals convergence ****

               *****************************************************
               *                     SUCCESS                       *
               *           SCF CONVERGED AFTER   1 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.5511e-17  Tolerance :   1.0000e-06
  Last MAX-Density change    ...    2.2204e-16  Tolerance :   1.0000e-05
  Last RMS-Density change    ...    4.7432e-17  Tolerance :   1.0000e-06
  Last DIIS Error            ...    2.5647e-16  Tolerance :   1.0000e-06

             **** THE GBW FILE WAS UPDATED (h.gbw) ****
             **** DENSITY FILE WAS UPDATED (h.scfp.tmp) ****
             **** ENERGY FILE WAS UPDATED (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 
--------------------------
 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
 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 0 sec 

Total time                  ....       0.478 sec
Sum of individual times     ....       0.445 sec  ( 93.0%)

Fock matrix formation       ....       0.285 sec  ( 59.7%)
Diagonalization             ....       0.000 sec  (  0.1%)
Density matrix formation    ....       0.000 sec  (  0.0%)
Population analysis         ....       0.001 sec  (  0.3%)
Initial guess               ....       0.150 sec  ( 31.4%)
Orbital Transformation      ....       0.000 sec  (  0.0%)
Orbital Orthonormalization  ....       0.000 sec  (  0.0%)
DIIS solution               ....       0.000 sec  (  0.0%)

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

MakePlots. Basename = [h]

---------------
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       ... 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       ... 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       ... 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       ... 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                                         ... h.gbw
Electron density file                           ... h.scfp.tmp
The origin for moment calculation is the CENTER OF MASS  = ( 0.000000,  0.000000  0.000000)

-------------
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         ...       13.986 sec (=   0.233 min)
GTO integral calculation        ...        4.723 sec (=   0.079 min)  33.8 %
SCF iterations                  ...        4.831 sec (=   0.081 min)  34.5 %
Orbital/Density plot generation ...        4.432 sec (=   0.074 min)  31.7 %
                             ****ORCA TERMINATED NORMALLY****
TOTAL RUN TIME: 0 days 0 hours 0 minutes 18 seconds 683 msec

Скрипт для отображения результатов в PyMOL:

with open('orbitals.pml', 'w+') as wf: wf.write('cmd.volume_ramp_new(\'ramp001\', [\ -0.015, 1.00, 0.00, 0.00, 0.00, \ -0.01, 1.00, 1.00, 0.00, 0.20, \ -0.005, 0.00, 0.00, 1.00, 0.00, \ 0.005, 0.00, 0.00, 1.00, 0.00, \ 0.01, 0.00, 1.00, 1.00, 0.20, \ 0.015, 0.00, 0.00, 1.00, 0.00, \ ])\n') wf.write('set_view (\ -0.627283990, 0.770159245, -0.115637451,\ 0.233201414, 0.044078715, -0.971426964,\ -0.743056655, -0.636328280, -0.207251951,\ 0.000000000, 0.000000000, -141.782043457,\ -15.499999046, -15.499999046, -15.499999046,\ 111.782043457, 171.782043457, -20.000000000 )\n') for n in range(0,4): for l in range(0,n): for m in range(0,l+1,1): name='%s-%s-%s' % (n,l,m) lname = name+'.cube' volname = name+'_vol' wf.write('load %s, %s\n' % (lname, name)) wf.write('volume %s, %s\n' % (volname, name)) wf.write('orbit %s, ramp001\n' % (volname))

Из выдачи ORCA были визуализированы только 2 и 3 уровни:

  • Первое изображение 2-0-0 2-1-0 2-1-1

  • Второе изображение 3-0-0 3-1-0 3-1-1 3-2-0 3-2-1 3-2-2

In [22]:
m = Image(filename='h_1_2.png')
n = Image(filename='h_3.png')
display(m, n)

Ссылки

In [ ]: