论坛 › 教程分享 › Material Studio之perl脚本与cif，blender之python建模等

tfsn20 发表于 2024-1-17 20:30:43

Material Studio之perl脚本与cif，blender之python建模等

本帖最后由 tfsn20 于 2024-2-3 17:25 编辑

😎目的：导出MS建模的原子和键坐标，用于blender中渲染（[更好看](https://zhuanlan.zhihu.com/p/28703836)）
# MS之cif（用作了解）
各个软件对导出的cif文件内容都有所不同，😀不同版本对cif文件的解析有所不同，cif也有自己的版本更新，软件版本为Materials Studio 2020 (20.1.0.2728)。
## 导出graphite的cif文件
😃打开MS，新建项目，导入MS自带的graphite.msi的文件，导出选项中选择保存类型为cif。
## 解析
可以参考[😀](https://journals.iucr.org/j/issues/2016/01/00/aj5269/index.html)和[😄](https://www.iucr.org/resources/cif/spec/version1.1/cifsyntax)对各项的解析。在解析之前，让我解释一下CIF的基本结构：

    Data Block（数据块）： CIF文件通常包含一个或多个数据块，每个数据块以 data_ 关键字开始，后跟一个唯一的标识符，如 data_my_structure。

    Data Items（数据项）： 数据块中包含数据项，它们描述了实验数据、晶体结构等信息。每个数据项都有一个唯一的标签。

    Loop Block（循环块）： 循环块用于以表格形式列出一系列相关的数据项。它以 loop_ 关键字开始，后跟一系列标签，表示在接下来的数据行中将列出这些标签对应的值。
😄用记事本打开该cif文件。
```
data_graphite
_audit_creation_date            2024-01-14
_audit_creation_method            'Materials Studio'
_symmetry_space_group_name_H-M    'P63/MMC'
_symmetry_Int_Tables_number       194
_symmetry_cell_setting            hexagonal
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
-y,x-y,z
-x+y,-x,z
-x,-y,z+1/2
y,-x+y,z+1/2
x-y,x,z+1/2
y,x,-z
x-y,-y,-z
-x,-x+y,-z
-y,-x,-z+1/2
-x+y,y,-z+1/2
x,x-y,-z+1/2
-x,-y,-z
y,-x+y,-z
x-y,x,-z
x,y,-z+1/2
-y,x-y,-z+1/2
-x+y,-x,-z+1/2
-y,-x,z
-x+y,y,z
x,x-y,z
y,x,z+1/2
x-y,-y,z+1/2
-x,-x+y,z+1/2
_cell_length_a                  2.4600
_cell_length_b                  2.4600
_cell_length_c                  6.8000
_cell_angle_alpha               90.0000
_cell_angle_beta                  90.0000
_cell_angle_gamma               120.0000
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_U_iso_or_equiv
_atom_site_adp_type
_atom_site_occupancy
C2   C   0.00000   0.00000   0.25000   0.00000Uiso   1.00
C4   C   0.33333   0.66667   0.25000   0.00000Uiso   1.00
loop_
_geom_bond_atom_site_label_1
_geom_bond_atom_site_label_2
_geom_bond_distance
_geom_bond_site_symmetry_2
_ccdc_geom_bond_type
C2   C4      1.420   .   A
C2   C4      1.420   1_445 A
C2   C4      1.420   1_545 A
C4   C2      1.420   1_665 A
C4   C2      1.420   1_565 A
```
😁第一个数据块
```
data_graphite
_audit_creation_date            2024-01-14
_audit_creation_method            'Materials Studio'
_symmetry_space_group_name_H-M    'P63/MMC'
_symmetry_Int_Tables_number       194
_symmetry_cell_setting            hexagonal
```
(https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Isymmetry_space_group_name_H-M.html)：晶体所属的对称群名称，被_space_group_name_H-M_alt取代。
(https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Isymmetry_Int_Tables_number.html)：晶体所属对称群名称的数字编号，被_space_group_IT_number取代。
(https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Isymmetry_cell_setting.html)：晶体所属的八大晶系名称，被 (https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Ispace_group_crystal_system.html)取代。
😊第二个循环块：
```
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
-y,x-y,z
-x+y,-x,z
-x,-y,z+1/2
y,-x+y,z+1/2
x-y,x,z+1/2
y,x,-z
x-y,-y,-z
-x,-x+y,-z
-y,-x,-z+1/2
-x+y,y,-z+1/2
x,x-y,-z+1/2
-x,-y,-z
y,-x+y,-z
x-y,x,-z
x,y,-z+1/2
-y,x-y,-z+1/2
-x+y,-x,-z+1/2
-y,-x,z
-x+y,y,z
x,x-y,z
y,x,z+1/2
x-y,-y,z+1/2
-x,-x+y,z+1/2
_cell_length_a                  2.4600
_cell_length_b                  2.4600
_cell_length_c                  6.8000
_cell_angle_alpha               90.0000
_cell_angle_beta                  90.0000
_cell_angle_gamma               120.0000
```
_symmetry_equiv_pos_as_xyz：晶体所属空间群的对称操作，可以利用这个算原子坐标，但我不会（主要涉及到群论，太专业了😅），被(https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Ispace_group_symop_operation_xyz.html)取代。
\_cell\_length\_\*，\_cell\_angle\_*：晶格常数，很有用😎。
🙂第三个是循环块：
8个数据key
```
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_U_iso_or_equiv
_atom_site_adp_type
_atom_site_occupancy
```
对应多组5列数据
```
C2   C   0.00000   0.00000   0.25000   0.00000Uiso   1.00
C4   C   0.33333   0.66667   0.25000   0.00000Uiso   1.00
```
🙃最后一个loop块同理。
但cif文件里并没有明显的原子的坐标信息（可以用群论解决）。
但在MS可以用perl脚本，导出所有原子坐标。
# MS之导出xsd文件中显示范围内的所有原子坐标和键的起始位置（单独原子键除外）
## 文件类型
有3DAtomistic，PerlScript，ModuleState，
```
#!perl

use strict;
use warnings;
use Getopt::Long;
use MaterialsScript qw(:all);


my $doc = $Documents{"graphite.xsd"};

# DisplayRange avert Cannot operate upon infinite SymmetrySystem (function/property "Atoms")
my $atoms = $doc->DisplayRange->Atoms;
my $atoms_pos = '';
my $attached_atoms_pos='';
foreach my $index (0..$#{$atoms}) {
    my $atom = $atoms->[$index];
    my $x = $atom->X;
    my $y = $atom->Y;
    my $z = $atom->Z;

    my $connectedAtoms = $atom->AttachedAtoms;
    foreach my $atom_ (@$connectedAtoms) {
      $attached_atoms_pos=$attached_atoms_pos.'['.$atom_->X.' '.$atom_->Y.' '.$atom_->Z.']';
    };
    $attached_atoms_pos=$attached_atoms_pos."\n";
    $atoms_pos = $atoms_pos.$atom->X.' '.$atom->Y.' '.$atom->Z."\n";
};
#print $atoms_pos;
#print $attached_atoms_pos;

my $report=Documents->New("pos.txt");
$report->Append($atoms_pos.$attached_atoms_pos);
$report->Close;

$doc->Discard;
```
😀pos.txt文件即可做下一步使用。
# blender之python建模
在代码中更改pos.txt文件路径
```
import re
# 打开文本文件
with open(r"D:\Documents\School\MS\defaultNameProject_Files\Documents\graphite CASTEP Energy\pos.txt", 'r') as file:
    # 读取每一行内容
    lines = file.readlines()
if lines[-1] == '\n':
    lines.pop()
#get txt rows
txt_rows_length=len(lines)
# 初始化空数组用于存储结果
atoms_pos = []
attached_atoms_pos=[]
# 遍历每一行内容
for i in range(0, int(txt_rows_length/2)):
    # 将每一行按空格分割，并将结果转换为浮点数
    line = lines
    values =
    # 将结果添加到数组
    atoms_pos.append(values)
for i in range(int(txt_rows_length/2), txt_rows_length):
    # 将每一行按空格分割，并将结果转换为浮点数
    line = lines
    # 使用正则表达式提取方括号内的内容
    matches = re.findall(r'\[([^\]]+)\]', line)
    values = [ for e in matches]
    # 将结果添加到数组
    attached_atoms_pos.append(values)
# 打印最终的数组
# print(atoms_pos)
# print(attached_atoms_pos)

bonds_link_digital_serial_number=[]
for i,a in enumerate(attached_atoms_pos) :
    for j,b in enumerate(a):
      for k,c in enumerate(atoms_pos):
            if c == b:
                bonds_link_digital_serial_number.append('-'.join(sorted(f'{i}-{k}'.split('-'))))

bonds_link_digital_serial_number=list(set(bonds_link_digital_serial_number))
print(bonds_link_digital_serial_number)
bonds_link_pos=[]

for i,a in enumerate(bonds_link_digital_serial_number):
    num_pos =
    bonds_link_pos.append(],atoms_pos]])

print(bonds_link_pos)
import bpy
import math
# 清空场景
bpy.ops.object.select_all(action='DESELECT')
bpy.ops.object.select_by_type(type='MESH')
bpy.ops.object.delete()

# 创建原子的坐标信息
atom_coordinates = atoms_pos

# 创建原子
for coord in atom_coordinates:
    bpy.ops.mesh.primitive_uv_sphere_add(radius=0.2, location=coord)

def cylinder_between(x1, y1, z1, x2, y2, z2, r):

dx = x2 - x1
dy = y2 - y1
dz = z2 - z1   
dist = math.sqrt(dx**2 + dy**2 + dz**2)

bpy.ops.mesh.primitive_cylinder_add(
      radius = r,
      depth = dist,
      location = (dx/2 + x1, dy/2 + y1, dz/2 + z1)   
)

phi = math.atan2(dy, dx)
theta = math.acos(dz/dist)

bpy.context.object.rotation_euler = theta
bpy.context.object.rotation_euler = phi


# 创建键
for coords in bonds_link_pos:
    cylinder_between(coords,coords,coords,coords,coords,coords,0.1)


# 切换到对象模式
bpy.ops.object.mode_set(mode='OBJECT')

```
粘贴代码运行（Atl+P快捷键运行）
😃完成建模

查看完整版本: Material Studio之perl脚本与cif，blender之python建模等