Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2017

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

Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами

Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol.

   1 import numpy
   2 import scipy.special
   3 import scipy.misc

   1 import npy2cube

   1 def w(n,l,m,d):
   2     
   3     x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
   4     
   5     r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
   6     theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
   7     phi = lambda x,y,z: numpy.arctan(y/x)
   8     
   9     a0 = 1.
  10     
  11     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)
  12     WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
  13     absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
  14 
  15     ### может тут чего то не хватает?  
  16     return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

   1 # определите шаг grid при заданном диапозоне от -d до d
   2 d=.....
   3 step= .....
   4 
   5 # Зададим цикл по перебору квантовых чисел
   6 
   7 for n in range(0,4):
   8     for l in range(0,n):
   9         for m in range(0,l+1,1):
  10             grid= .... 
  11             name='%s-%s-%s' % (n,l,m)
  12             # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
  13             npy2cube.npy2cube(grid,(-d,....),(step,.....),name+'.cube')

import py3Dmol
view = py3Dmol.view()
alpha = open('1-0-0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.04, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()

   1 ### Откуда эти цифры? 
   2 pml='''
   3 ### cut below here and paste into script ###
   4 cmd.volume_ramp_new('ramp007', [\
   5       0.005, 0.00, 0.00, 1.00, 0.00, \
   6       0.01,  0.00, 1.00, 1.00, 0.20, \
   7       0.015, 0.00, 0.00, 1.00, 0.00, \
   8     ])
   9 ### cut above here and paste into script ###
  10 '''
  11 
  12 for n in range(0,4):
  13     for l in range(0,n):
  14         for m in range(0,l+1,1):
  15             name='%s-%s-%s' % (n,l,m)
  16             # Загрузить cube файл
  17             # Отобразить электронную плотность (hint: volume)
  18             # Покрасить ее разумным образом 
  19 
  20 # запишите переменную в файл
  21 v=open(.....

   1 from IPython.display import display,Image
   2 display(Image=(file))

g = '''
         0 1
         C     0.000000     0.000000     0.000000
'''

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

2017/8/task4 (последним исправлял пользователь golovin 2023-03-10 11:30:44)