In [ ]:
В ходе выполнения данного задания столкнулся с проблемой, как оказалось PyMol не может адекватно построить volume вокруг грида обычной волновой функции,
поэтому пришлось прибегнуть к небольшой хитрости и строить объёмы вокруг абсолютной волновой функции.
In [ ]:
import npy2cube
In [ ]:
import numpy
import scipy.special
import scipy.misc
In [ ]:
def w(n,l,m,d):
    
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    
    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.
    
    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))

    ### может тут чего то не хватает?  
    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)#, WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
In [ ]:
# определите шаг grid при заданном диапозоне от -d до d
d=20
step=(0.5, 0.5, 0.5)

# Зададим цикл по перебору квантовых чисел

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,name+'.cube')
In [86]:
from IPython.display import Image
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]


import pymol
pymol.finish_launching()
from pymol import cmd,stored



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)
            cmd.load(name+'.cube', 'orbital_map')
            cmd.volume('orbital_volume', 'orbital_map')
            cmd.zoom('orbital_map', buffer=5)
            cmd.orient('orbital_map')
            cmd.set('fog', 0)
            cmd.turn('y', 90)
            cmd.turn('z', 45)
            cmd.turn('x', 20)
            cmd.png(name+'.png')
            cmd.reinitialize()
In [87]:
Image(filename='1-0-0.png')
Out[87]:
No description has been provided for this image
In [88]:
Image(filename='2-0-0.png')
Out[88]:
No description has been provided for this image
In [89]:
Image(filename='2-1-0.png')
Out[89]:
No description has been provided for this image
In [90]:
Image(filename='2-1-1.png')
Out[90]:
No description has been provided for this image
In [91]:
Image(filename='3-0-0.png')
Out[91]:
No description has been provided for this image
In [92]:
Image(filename='3-1-0.png')
Out[92]:
No description has been provided for this image
In [93]:
Image(filename='3-1-1.png')
Out[93]:
No description has been provided for this image
In [94]:
Image(filename='3-2-0.png')
Out[94]:
No description has been provided for this image
In [95]:
Image(filename='3-2-1.png')
Out[95]:
No description has been provided for this image
In [96]:
Image(filename='3-2-2.png')
Out[96]:
No description has been provided for this image

Теперь выполним второе задание и рассчитаем орбитали с помощью программы Psi4:

In [ ]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [ ]:
g = '''
         0 1
         C     0.000000     0.000000     0.000000
'''
In [ ]:
m = psi4.geometry(g)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
ener=psi4.energy('scf/cc-pvtz', molecule = m )
#Расчёт орбиталей в Psi4:  информация  об api https://psicode.org/psi4manual/1.4.0/api/psi4.core.cubeproperties
In [ ]:
molecule = psi4.geometry(g)
psi4.set_options({
    'basis': 'cc-pVDZ',           
    'scf_type': 'df',              
    'freeze_core': True,           
    'cubeprop_tasks': ['orbitals'], 
    'cubeprop_orbitals': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 
    'cubic_grid_spacing': [0.2, 0.2, 0.2], 
    'cubic_grid_overage': [4.0, 4.0, 4.0]  
})

scf_energy, wavefunction = psi4.energy('scf', return_wfn=True, molecule=molecule)
psi4.cubeprop(wavefunction)
In [1]:
from IPython.display import Image
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]


import pymol
pymol.finish_launching()
from pymol import cmd,stored
In [109]:
for n in range(0,11):
    name=str(n)
    cmd.load(name+'.cube', 'orbital_map')
    cmd.enable('orbital_map')
    cmd.hide('extent')
    cmd.show('dots')
    cmd.set('fog', 0)
    cmd.turn('x', 90)
    cmd.turn('z', 55)
    cmd.turn('y', 35)
    cmd.png(name+'_orbital.png')
    cmd.reinitialize()

К сожалению, если брать не абсолютное значение волновой функции, нормально построить объём вокруг имеющегося грида паймол не может. Поэтому визуализируем посчитанные psi4 орбитали с помощью dots.

In [110]:
Image(filename='1_orbital.png')
Out[110]:
No description has been provided for this image
In [111]:
Image(filename='2_orbital.png')
Out[111]:
No description has been provided for this image
In [112]:
Image(filename='3_orbital.png')
Out[112]:
No description has been provided for this image
In [113]:
Image(filename='4_orbital.png')
Out[113]:
No description has been provided for this image
In [114]:
Image(filename='5_orbital.png')
Out[114]:
No description has been provided for this image
In [115]:
Image(filename='6_orbital.png')
Out[115]:
No description has been provided for this image
In [116]:
Image(filename='7_orbital.png')
Out[116]:
No description has been provided for this image
In [117]:
Image(filename='8_orbital.png')
Out[117]:
No description has been provided for this image
In [118]:
Image(filename='9_orbital.png')
Out[118]:
No description has been provided for this image
In [119]:
Image(filename='10_orbital.png')
Out[119]:
No description has been provided for this image
In [ ]: