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

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

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

😎目的：导出MS建模的原子和键坐标，用于blender中渲染（更好看）

MS之cif（用作了解）

各个软件对导出的cif文件内容都有所不同，😀不同版本对cif文件的解析有所不同，cif也有自己的版本更新，软件版本为Materials Studio 2020 (20.1.0.2728)。

导出graphite的cif文件

😃打开MS，新建项目，导入MS自带的graphite.msi的文件，导出选项中选择保存类型为cif。
图片.png

解析

可以参考😀和😄对各项的解析。在解析之前，让我解释一下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.00000  Uiso   1.00
C4     C     0.33333   0.66667   0.25000   0.00000  Uiso   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

_symmetry_space_group_name_H-M：晶体所属的对称群名称，被_space_group_name_H-M_alt取代。
_symmetry_Int_Tables_number：晶体所属对称群名称的数字编号，被_space_group_IT_number取代。
_symmetry_cell_setting：晶体所属的八大晶系名称，被 _space_group_crystal_system取代。
😊第二个循环块：

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：晶体所属空间群的对称操作，可以利用这个算原子坐标，但我不会（主要涉及到群论，太专业了😅），被_space_group_symop_operation_xyz取代。
_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.00000  Uiso   1.00
C4     C     0.33333   0.66667   0.25000   0.00000  Uiso   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;

图片.png

😀pos.txt文件即可做下一步使用。

blender之python建模

图片.png
在代码中更改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[i]
    values = [float(x) for x in line.split()]
    # 将结果添加到数组
    atoms_pos.append(values)
for i in range(int(txt_rows_length/2), txt_rows_length):
    # 将每一行按空格分割，并将结果转换为浮点数
    line = lines[i]
    # 使用正则表达式提取方括号内的内容
    matches = re.findall(r'\[([^\]]+)\]', line)
    values = [[float(x) for x in e.split()] 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 = [int(num) for num in a.split('-')]
    bonds_link_pos.append([atoms_pos[num_pos[0]],atoms_pos[num_pos[1]]])

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[1] = theta 
  bpy.context.object.rotation_euler[2] = phi 

# 创建键
for coords in bonds_link_pos:
    cylinder_between(coords[0][0],coords[0][1],coords[0][2],coords[1][0],coords[1][1],coords[1][2],0.1)

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

粘贴代码运行（Atl+P快捷键运行）
😃完成建模
图片.png