计算化学公社

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

---

作者
Author: 鬼隐    时间: 2022-3-7 17:00
标题: Blender + Python:用少量有效数据绘制势能面示意图的方法

前几天有需求要绘制一种势能面的示意图，类似教科书上标出一阶鞍点、 局域极小点那种示意图。

这种图正规绘制需要大量的单点计算，并用软件描面画成。但是实际上，我无法计算出如此多的单点来绘制一张图，毕竟在一般的计算有机工作中，单是定位过渡态就够普通鼠标侠喝一壶了，更别说选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.

• 首先下载安装Blender，然后打开软件，可以更改下语言为中文。
• 删除初始产生的对象，最多只留下光源和摄像机
• 点击 添加 >> 网格 >> 平面
  
  点这里
• 点击上图中的 物体模式改为编辑模式
  随后左边多出一竖列按钮，鼠标移到上面会显示名称，点击环切
  同时编辑模式字样下方出现切割次数，建议改为10或者8,9
• 鼠标移到添加的平面上，会出现一根黄线，点击，垂直各切一次
  
  点这里
  
  点这里
• 鼠标左键长按该列第一个按钮，选中刷选
  
  点这里
• 平面上切出的格子点上按住鼠标移动，选中
• 按 G便可以自由移动选中的点，建议再按Z限制其只改变Z坐标,按住鼠标中键可改变视角。
  （当然，我们可以使用调整，直接捏出想要的形状，但是我们目的不是在Blender中捏出，只是利用其捏出基本形貌，后利用坐标在常用的科研作图软件中插值制作出符合审美的示意图，所以我建议使用刷选并只改变Z值，你也可以衰减编辑, 但我觉得操作单个点舒服，而且只变Z值的话，xy坐标均匀，后续代码拟合曲面效果会好）
  这些点应该严格按照已有的数据绘制出相对高低，峰、谷、鞍面合理处理位置。
• 使用鼠标拖动，竖直改变对象的形貌，大致捏出基本形貌，该低的地方低，该高的地方高，哪儿是鞍点那儿就捏成鞍部。
  比如我可以随手捏出这么一个ts连接两个minimum的示意图，看起来比较丑陋
  
  

捏平面
导出为xyz坐标

• 此时保存一下自己捏的势能草面，然后进入物体模式，鼠标拖拉选中势能草面
• 文件 >> 导出 >> Waveforont(.obj)
  勾中仅导出选中的物体
  几何数据只选三角面
  很幸运，.obj是文本可以进行文本解析，也可以用windows自带的3D查看器打开
  
  3D查看器
• 提取坐标
  打开git bash然后输入

  1. grep "v " surface.obj |awk '{print $2,$4,$3}' |sort > surface.xyz

  复制代码

  
  需要注意的是，导出.obj时坐标可以变换，在不改动设置的情况下，导出的坐标是y值和z值放反了，于是用awk 换下顺序，并按照x值大小排序。具体变换可以自己操作
  
  

使用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. 
   
7. with open("surface.xyz","r") as f:
   
8.     s = f.read()
   
9.    
   
10. xyz = s.splitlines()
    
11. 
    
12. x = []
    
13. y = []
    
14. z = []
    
15. for inf in range(len(xyz)):
    
16.     x.append(float(xyz[inf].split()[0]))
    
17.     y.append(float(xyz[inf].split()[1]))
    
18.     z.append(float(xyz[inf].split()[2]))
    
19. print(max(x),min(x),max(y),min(y),max(z),min(z))
    
20. 
    
21. dim = int(np.sqrt(len(xyz)))
    
22. X = np.asarray(x).reshape([dim, dim])
    
23. Y = np.asarray(y).reshape([dim, dim])
    
24. Z = np.asarray(z).reshape([dim, dim])
    
25. 
    
26. #fig = plt.figure()  
    
27. #ax = plt.axes(projection='3d')
    
28. #plt.axis('off')
    
29. #ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1, cmap = plt.get_cmap('gnuplot'))
    
30. #plot.show()
    
31. 
    
32. 
    
33. x_0 = np.linspace(min(x),max(x),200)
    
34. y_0 = np.linspace(min(y),max(y),200)
    
35. X_0, Y_0 = np.meshgrid(x_0, y_0,indexing='ij')
    
36. f = interpolate.Rbf(X, Y, Z, function='multiquadric',smooth=0.3)
    
37. Z_1 = f(X_0, Y_0)
    
38. Z_1.shape
    
39. 
    
40. 
    
41. fig = plt.figure()  
    
42. ax = plt.axes(projection='3d')
    
43. plt.axis('off')
    
44. 
    
45. ax.plot_surface(X_0, Y_0, Z_1, rstride = 1, cstride = 1, color="#cccccc")
    
46. 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


一个 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: 94)

上传 Uploaded

点击下载
Click to download

---

作者
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