Практикум 4¶
In [1]:
import numpy
import scipy.special
import scipy.misc
In [6]:
import npy2cube
Определим функцию w: (квантовые числа, объем) -> (сетка значений волновой функции)
In [4]:
from math import factorial
def w(n,l,m,d):
x,y,z = numpy.mgrid[-d:d:100j,-d:d:100j,-d:d:100j]
#переход в полярные координаты: выражение r, theta, phi (полярные) через x, y, z (декартовы)
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.
# запишем волновую фнкцию в сферических координатах (слайд 12 в лекции)
# обозначим L -- обобщённый полином Лагерра степени n-l-1
# обозначим sph -- сферическая гармоника (зависимость волновой функции от направления)
# R = e^(-ro/2) * ro^l * L
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)
# coeff = (2/(n*a0))**1.5 * ((n-l-1)!/(2n(n+l)!))^0.5
coeff = lambda n,l: ((factorial(n-l-1)/(2*n*factorial(n+l)))*(2/n/a0)**3) ** 0.5
# WF = R * sph * coeff
WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta) * coeff(n, l)
# квадрат модуля волновой функции (плотность вероятности нахождения электрона в точке),
# в этой функции не нужен
##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 [7]:
# определите шаг grid при заданном диапозоне от -d до d
d = 30
step= (2 * d + 1) / 100
# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
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)
# для сохранения нужно задать координаты старта grid и шаг по каждому направлению
npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
In [2]:
from IPython.display import HTML
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
Создадим скрипт PyMOL, отображающий полученные орбитали.
In [8]:
with open('volume_color.pml', 'w') as out:
# ramp -- отображение (значение WF) -> (цвет в PyMOL)
# оно задано по точкам в виде пятерок {значение, R, G, B, alpha}
# (здесь отображение определено в 6 точках, в т. ч. отрицательных)
out.write('cmd.volume_ramp_new(\'ramp007\', [\
-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')
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)
cubename = name+'.cube'
volname = name+'_vol'
out.write('load {0}, {1}\n'.format(cubename, name))
out.write('volume {0}, {1}\n'.format(volname, name))
out.write('volume_color {0}, ramp007\n'.format(volname))
In [16]:
cmd.load("volume_color.pml")
cmd.bg_color('white')
In [3]:
cmd.zoom()
Сделаем картинки полученных орбиталей в PyMOL (вручную) и отобразим их
In [15]:
import matplotlib.pyplot as plt
import glob
from os.path import basename, splitext
f, axarr = plt.subplots(5, 2, figsize=(20,20))
curr_row = 0
for n, f in enumerate(sorted(glob.glob('pr4/py/*.png'))):
# fetch the url as a file type object, then read the image
a = plt.imread(f)
col = n % 2
# plot on relevant subplot
axarr[curr_row, col].imshow(a)
axarr[curr_row, col].axis('off')
axarr[curr_row, col].set_title(basename(f))
if col == 1:
# we have finished the current row, so increment row counter
curr_row += 1
Создадим входной файл для Orca:
In [22]:
with open('h.inp', 'w') as f:
f.write('''! 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
*''')
Запустим на kodomo: 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 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 = 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.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: 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.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
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.000 sec
Total time needed ... 0.001 sec
---------------------
INITIAL GUESS: MOREAD
---------------------
Guess MOs are being read from file: h.ges
Input Geometry matches current geometry (good)
Input basis set matches current basis set (good)
MOs were renormalized
MOs were reorthogonalized (Cholesky)
------------------
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 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 ... 1.6653e-16 Tolerance : 1.0000e-06
Last MAX-Density change ... 2.2204e-16 Tolerance : 1.0000e-05
Last RMS-Density change ... 4.6693e-17 Tolerance : 1.0000e-06
Last DIIS Error ... 7.3388e-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
-----------------------------------------------
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.441675046686;
(*BG)[ 1][ 2].a[ 1]= 0.80000000; (*BG)[ 1][ 2].d[ 1]=-0.637945909953;
(*BG)[ 1][ 2].a[ 2]= 0.80000000; (*BG)[ 1][ 2].d[ 2]= 0.749010191077;
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.409543739990
2 0.800000000000 -0.591536143673
3 0.800000000000 0.694520637392
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.430 sec
Sum of individual times .... 0.412 sec ( 95.8%)
Fock matrix formation .... 0.373 sec ( 86.7%)
Diagonalization .... 0.001 sec ( 0.1%)
Density matrix formation .... 0.000 sec ( 0.0%)
Population analysis .... 0.020 sec ( 4.8%)
Initial guess .... 0.018 sec ( 4.1%)
Orbital Transformation .... 0.000 sec ( 0.0%)
Orbital Orthonormalization .... 0.018 sec ( 4.1%)
DIIS solution .... 0.000 sec ( 0.0%)
------------------------- --------------------
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 ... 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
-------------
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 ... 1.156 sec (= 0.019 min)
GTO integral calculation ... 0.329 sec (= 0.005 min) 28.5 %
SCF iterations ... 0.468 sec (= 0.008 min) 40.5 %
Orbital/Density plot generation ... 0.359 sec (= 0.006 min) 31.0 %
****ORCA TERMINATED NORMALLY****
TOTAL RUN TIME: 0 days 0 hours 0 minutes 1 seconds 373 msecАналогично визуализируем орбитали:
In [10]:
with open('volume_color_orca.pml', 'w') as out:
# ramp -- отображение (значение WF) -> (цвет в PyMOL)
# оно задано по точкам в виде пятерок {значение, R, G, B, alpha}
# (здесь отображение определено в 6 точках, в т. ч. отрицательных)
out.write('cmd.volume_ramp_new(\'ramp007\', [\
-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')
for n in sorted(glob.glob('pr4/orca/*.cube')):
name=splitext(basename(n))[0]
cubename = name+'.cube'
volname = name+'_vol'
out.write('load {0}, {1}\n'.format(cubename, name))
out.write('volume {0}, {1}\n'.format(volname, name))
out.write('volume_color {0}, ramp007\n'.format(volname))
In [11]:
cmd.load("volume_color_orca.pml")
cmd.bg_color('white')
In [ ]:
cmd.zoom()
In [13]:
import matplotlib.pyplot as plt
import glob
from os.path import basename, splitext
f, axarr = plt.subplots(2,2, figsize=(20,20))
curr_row = 0
for n, f in enumerate(sorted(glob.glob('pr4/orca/*.png'))):
# fetch the url as a file type object, then read the image
a = plt.imread(f)
col = n % 2
# plot on relevant subplot
axarr[curr_row, col].imshow(a)
axarr[curr_row, col].axis('off')
axarr[curr_row, col].set_title(basename(f))
if col == 1:
# we have finished the current row, so increment row counter
curr_row += 1
Результат работы orca -- четыре орбитали, одна аналогична орбитали с квантовыми числами 1-0-0, три другие -- орбитали 2-1-0. При этом разрешение моделирования довольно низкое: в орбиталях 2-1-0 orca нельзя (не только на этих картинках с не очень высоким разрешением, но и в pymol) увидеть, что между "половинками" есть просвет.