上一主题 下一主题
ScriptCat,新一代的脚本管理器脚本站,与全世界分享你的用户脚本油猴脚本开发指南教程目录
返回列表 发新帖
楼主: tfsn20 - 

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

[复制链接]
  • TA的每日心情
    慵懒
    11 小时前
  • 签到天数: 765 天

    [LV.10]以坛为家III

    46

    主题

    197

    回帖

    849

    积分

    荣誉开发者

    积分
    849

    荣誉开发者油中2周年生态建设者

    发表于 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
    图片.png
    图片.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

    该用户从未签到

    0

    主题

    1

    回帖

    2

    积分

    助理工程师

    积分
    2
    发表于 2024-6-30 15:50:01 | 显示全部楼层
    这个创建键有点问题,难道是我模型太复杂?有的键特别长,有的穿模,有的堆在一起。。。另外颜色也是个问题,ms中颜色是一半对一半的,不知道有没有什么好的处理方式
    回复

    使用道具 举报

    发表回复

    本版积分规则

    快速回复 返回顶部 返回列表