计算化学公社

标题: Blender + Python:用少量有效数据绘制势能面示意图的方法 [打印本页]

作者
Author:
鬼隐    时间: 2022-3-7 17:00
标题: Blender + Python:用少量有效数据绘制势能面示意图的方法
本帖最后由 鬼隐 于 2022-3-8 10:31 编辑

前几天有需求要绘制一种势能面的示意图,类似教科书上标出一阶鞍点、 局域极小点那种示意图。
这种图正规绘制需要大量的单点计算,并用软件描面画成。但是实际上,我无法计算出如此多的单点来绘制一张图,毕竟在一般的计算有机工作中,单是定位过渡态就够普通鼠标侠喝一壶了,更别说选CV扫描并绘图了。时间更不允许,当我需要画图时,往往意味着科研需要的数据已经满足需求了,我再要求提供给我更多数据也不大可能得到。
因此只能利用现用的计算数据合理规划着绘制,也就是说我只能依据目前算出的几个结构的数据,或者IRC上的点进行绘制。
经朋友启发,组织了一个简易的流程,来画这种简易的示意图。
使用Blender绘制草稿
Blender是一款开源免费的建模软件
Blender is a free and open-source 3D computer graphics software toolset used for creating animated films, visual effects, art, 3D printed models, motion graphics, interactive 3D applications, virtual reality, and computer games. Blender's features include 3D modelling, UV unwrapping, texturing, raster graphics editing, rigging and skinning, fluid and smoke simulation, particle simulation, soft body simulation, sculpting, animating, match moving, rendering, motion graphics, video editing, and compositing.

捏平面
导出为xyz坐标使用Python或者其他软件处理坐标
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from scipy.interpolate import make_interp_spline
  4. from sklearn.kernel_ridge import KernelRidge
  5. from scipy import interpolate

  6. with open("surface.xyz","r") as f:
  7.     s = f.read()
  8.    
  9. xyz = s.splitlines()

  10. x = []
  11. y = []
  12. z = []
  13. for inf in range(len(xyz)):
  14.     x.append(float(xyz[inf].split()[0]))
  15.     y.append(float(xyz[inf].split()[1]))
  16.     z.append(float(xyz[inf].split()[2]))
  17. print(max(x),min(x),max(y),min(y),max(z),min(z))

  18. dim = int(np.sqrt(len(xyz)))
  19. X = np.asarray(x).reshape([dim, dim])
  20. Y = np.asarray(y).reshape([dim, dim])
  21. Z = np.asarray(z).reshape([dim, dim])

  22. #fig = plt.figure()  
  23. #ax = plt.axes(projection='3d')
  24. #plt.axis('off')
  25. #ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1, cmap = plt.get_cmap('gnuplot'))
  26. #plot.show()


  27. x_0 = np.linspace(min(x),max(x),200)
  28. y_0 = np.linspace(min(y),max(y),200)
  29. X_0, Y_0 = np.meshgrid(x_0, y_0,indexing='ij')
  30. f = interpolate.Rbf(X, Y, Z, function='multiquadric',smooth=0.3)
  31. Z_1 = f(X_0, Y_0)
  32. Z_1.shape


  33. fig = plt.figure()  
  34. ax = plt.axes(projection='3d')
  35. plt.axis('off')

  36. ax.plot_surface(X_0, Y_0, Z_1, rstride = 1, cstride = 1, color="#cccccc")
  37. plt.show()</code></pre>
复制代码

在被注释的 plt.show()那里会展示最初捏出的草图形貌,不过在Blender里就能看到它就不需要matplotlib渲染了
代码默认生成200x200的图,使用Rbf来拟合曲面。一切都可以改。代码就在上面
也可以保存坐标使用MatLab,Origin等专业软件进一步处理。
结果展示
刚才随意搓的示意图感觉还凑合,还是有些别扭,可能需要再调一调。
刚才的示意图
之前画的几个图
首先是捏的草图
最后调cmap = plt.get_cmap('gnuplot')着色的示意图

本文使用Zhihu on VScode 发布至知乎,随后复制网页到论坛编辑器,发现大部分格式都保留了,sob老师,坛子现在可以用markdown发帖吗?
作者
Author:
万里云    时间: 2022-3-7 17:31
这个能用高斯函数去近似吧?
作者
Author:
鬼隐    时间: 2022-3-7 17:37
万里云 发表于 2022-3-7 17:31
这个能用高斯函数去近似吧?

可以的,我想手捏可能容易调些,也比较有趣。
作者
Author:
喵星大佬    时间: 2022-3-7 17:45
http://bbs.keinsci.com/thread-19519-1-1.html
某神贴
作者
Author:
王二葛    时间: 2022-3-7 23:07
本帖最后由 王二葛 于 2022-3-7 23:09 编辑

一个 Mathematica 的小例子,供参考
  1. size = 5;
  2. data = Flatten[Table[{{x/size, y/size}, RandomReal[]}, {x, 0, size}, {y, 0, size}], 1];
  3. intp = Interpolation[data, Method -> "Spline"];
  4. Plot3D[intp[x, y], {x, 0, 1}, {y, 0, 1},
  5. ColorFunction -> "SunsetColors",
  6. ColorFunction -> Function[{x, y, z}, Hue[1 - z]], Boxed -> False,
  7. Axes -> False, Mesh -> 50, MeshStyle -> Gray]
复制代码

运行结果:
(, 下载次数 Times of downloads: 84)



作者
Author:
鬼隐    时间: 2022-3-7 23:29
万里云 发表于 2022-3-7 17:31
这个能用高斯函数去近似吧?

python中interpolate.Rbf()可以指定 function="gaussian", 代码处理的部分,方法多样,也可以griddata,或者interp2d。


我觉得这里比较方便的是,blender搓草图然后提取数据,可以用鼠标全部完成。
得到的xyz可以再单独操作具体数值,以符合计算结果的数据。
所以我建议环切,切大约8刀,得到10x10的格子,这样需要手动调整100个Z值的子集也轻松。






欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3