计算化学公社

 找回密码 Forget password
 注册 Register
Views: 1853|回复 Reply: 5

[数据作图] Blender + Python:用少量有效数据绘制势能面示意图的方法

[复制链接 Copy URL]

20

帖子

0

威望

1046

eV
积分
1066

Level 4 (黑子)

发表于 Post on 2022-3-7 17:00:59 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 鬼隐 于 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.
  • 首先下载安装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. 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发帖吗?

评分 Rate

参与人数
Participants 6
eV +23 收起 理由
Reason
muuu2333 + 4
vantasyoscar + 4 とてもいい!
biogon + 5 GJ!
hdhxx123 + 5 GJ!
zjxitcc + 2 好物!
hebrewsnabla + 3

查看全部评分 View all ratings

347

帖子

4

威望

2583

eV
积分
3010

Level 5 (御坂)

发表于 Post on 2022-3-7 17:31:13 | 显示全部楼层 Show all
这个能用高斯函数去近似吧?

评分 Rate

参与人数
Participants 1
eV +2 收起 理由
Reason
卡开发发 + 2 让我想起了matlab有高斯这个梗23333

查看全部评分 View all ratings

20

帖子

0

威望

1046

eV
积分
1066

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-3-7 17:37:23 | 显示全部楼层 Show all
万里云 发表于 2022-3-7 17:31
这个能用高斯函数去近似吧?

可以的,我想手捏可能容易调些,也比较有趣。

1464

帖子

1

威望

2599

eV
积分
4083

Level 6 (一方通行)

喵星人

发表于 Post on 2022-3-7 17:45:27 | 显示全部楼层 Show all

评分 Rate

参与人数
Participants 4
eV +11 收起 理由
Reason
qwoop + 3 承包本日笑点
muuu2333 + 4
鬼隐 + 3 乐死了,这下省了200块。
卡开发发 + 1 你太可爱

查看全部评分 View all ratings

57

帖子

1

威望

2531

eV
积分
2608

Level 5 (御坂)

发表于 Post on 2022-3-7 23:07:06 | 显示全部楼层 Show all
本帖最后由 王二葛 于 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]
复制代码

运行结果:

example_execution

example_execution



评分 Rate

参与人数
Participants 1
eV +2 收起 理由
Reason
鬼隐 + 2 wolfram引擎也免费了,很适合绘图

查看全部评分 View all ratings

十八介姑娘一蕾花呀,白白介牙齿、红红介嘴唇,得人惜

20

帖子

0

威望

1046

eV
积分
1066

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-3-7 23:29:53 | 显示全部楼层 Show all
万里云 发表于 2022-3-7 17:31
这个能用高斯函数去近似吧?

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


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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 02:59 , Processed in 0.352292 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list