对Lammps仿真结果分析的Python代码

发表于 更新于 分类于

数据读取函数

输出指定时间步的所有信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def readfiles(filename):
    with open(filename, 'r') as file:
        lines = file.readlines()
    return lines

def split_lammpstrj(filename, index): 
    lines = readfiles(filename)
    # result = []
    start_index = 0
    end_index = 0
    for i in range(len(lines)):
        if lines[i].strip() == "ITEM: TIMESTEP":
            #print(i)
            start_index = i + 1  # 下一行是单元的起始行
            end_index = start_index
            #print(start_index)
            
            if lines[start_index].strip() == str(index):
                #print(lines[start_index].strip())
                while 1:
                    #print(lines[end_index], "99")
                    if lines[end_index].strip() == "ITEM: TIMESTEP":
                        end_index = end_index - 1
                        break
                    elif end_index+1 == len(lines):
                        end_index = end_index
                        break
                    else:
                        
                        end_index += 1
                
                
                result = lines[start_index:end_index+1]
                break
    return result

获取模拟的设置信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def getSimulationSet(filename):

    target1 = "ITEM: NUMBER OF ATOMS"  # 目标字符串1
    target2 = "ITEM: BOX BOUNDS pp pp pp"  # 目标字符串2

    NumberAtoms = None  # 原子数
    bounds = []  # 边界

    with open(filename, 'r') as file:
        lines = file.readlines()
        found1 = False
        found2 = False
    
        for i, line in enumerate(lines):
            if target1 in line:
                NumberAtoms = int(lines[i+1].strip())
                found1 = True
            elif target2 in line:
                bounds.append(lines[i+1].strip())
                bounds.append(lines[i+2].strip())
                bounds.append(lines[i+3].strip())
                found2 = True
    
            if found1 and found2:
                break

    # 创建一个空的列表来存储结果
    bounds_list = []
    
    # 遍历每一行
    for line in bounds:
        # 使用空格进行分割，得到数字字符串列表
        number_str_list = line.split()
        
        # 将数字字符串列表转换为数字列表
        number_list = [float(num_str) for num_str in number_str_list]
        
        # 将每一行的数字列表添加到结果列表
        bounds_list.append(number_list)

    return NumberAtoms, np.array(bounds_list)

获取指定时间步的原子id、类型、坐标和速度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def getAtomsInformation(filename, index):
    data = split_lammpstrj(filename, index)
    
    output = []
    found_match = False
    
    for line in data:
        if "ITEM: ATOMS" in line:
            found_match = True
        elif found_match:
            output.append(line)
    
    string_list = output
    
    float_array_list = []

    for string in string_list:
        split_string = string.rstrip().split(' ')
        float_array = np.array(split_string, dtype=float)
        float_array_list.append(float_array)
    
    float_array_list = np.array(float_array_list)
    return float_array_list

获取所有时间步并存放到列表中

1
2
3
4
def getAllTimesteps(filename):
    with open(filename, 'r') as file:
        string_list = file.readlines()
    float_list = []

计算指定时间步、指定区域内坐标的原子数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def calculateTheAtomsNumbers(filename, index, ranges):
    x_min = ranges[0][0]
    x_max = ranges[0][1]
    y_min = ranges[1][0]
    y_max = ranges[1][1]
    z_min = ranges[2][0]
    z_max = ranges[2][1]
    
    N = 0
    ################################
    data = getAtomsInformation(filename, index)
    
    # Extract x, y, z coordinates from the data
    x2 = data[:, 2]
    y2 = data[:, 3]
    z2 = data[:, 4]
    
    for i in range(len(x2)):
        if x_min <= x2[i] < x_max:      
            if y_min <= y2[i] <= y_max: 
                if z_min <= z2[i] <= z_max:
    
                    N+=1
    return N
```    

## 绘制结果

### 根据返回的值绘制模拟边界(长方体边界)
```python
def plotBounds(filename):
    numberAtoms, bounds =  getSimulationSet(filename)
    
    x1 = bounds[0,0]
    x2 = bounds[0,1]
    y1 = bounds[1,0]
    y2 = bounds[1,1]
    z1 = bounds[2,0]
    z2 = bounds[2,1]



    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    # 计算长方体的8个顶点坐标
    x = [x1, x1, x1, x1, x2, x2, x2, x2]
    y = [y1, y1, y2, y2, y1, y1, y2, y2]
    z = [z1, z2, z2, z1, z1, z2, z2, z1]
    
    # 定义连接顶点的线段
    edges = [
        [0, 1], [1, 2], [2, 3], [3, 0],  # 底部四边
        [0, 4], [1, 5], [2, 6], [3, 7],  # 底部到顶部的四条垂直边
        [4, 5], [5, 6], [6, 7], [7, 4]   # 顶部四边
    ]
    
    # 绘制长方体边框线
    for edge in edges:
        x_edge = [x[edge[0]], x[edge[1]]]
        y_edge = [y[edge[0]], y[edge[1]]]
        z_edge = [z[edge[0]], z[edge[1]]]
        ax.plot(x_edge, y_edge, zs=z_edge, color='black')
    
    # 设置坐标轴范围
    ax.set_xlim([min(x), max(x)])
    ax.set_ylim([min(y), max(y)])
    ax.set_zlim([min(z), max(z)])
    
    # 关闭坐标轴网格线
    ax.grid(False)
    
    # 设置按照实际数值显示比例
    ax.axis('equal')
    
    # 隐藏坐标轴
    ax.axis('off')
    
    # 显示图形
    plt.show()

绘制指定时间步的原子分布

性能差

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def plotAtomsSites(filename,index, radius):
    
    numberAtoms, bounds =  getSimulationSet(filename)
    
    x1 = bounds[0,0]
    x2 = bounds[0,1]
    y1 = bounds[1,0]
    y2 = bounds[1,1]
    z1 = bounds[2,0]
    z2 = bounds[2,1]

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    # 计算长方体的8个顶点坐标
    X = [x1, x1, x1, x1, x2, x2, x2, x2]
    Y = [y1, y1, y2, y2, y1, y1, y2, y2]
    Z = [z1, z2, z2, z1, z1, z2, z2, z1]
    
    # 定义连接顶点的线段
    edges = [
        [0, 1], [1, 2], [2, 3], [3, 0],  # 底部四边
        [0, 4], [1, 5], [2, 6], [3, 7],  # 底部到顶部的四条垂直边
        [4, 5], [5, 6], [6, 7], [7, 4]   # 顶部四边
    ]
    
    # 绘制长方体边框线
    for edge in edges:
        x_edge = [X[edge[0]], X[edge[1]]]
        y_edge = [Y[edge[0]], Y[edge[1]]]
        z_edge = [Z[edge[0]], Z[edge[1]]]
        ax.plot(x_edge, y_edge, zs=z_edge, color='black')
    

    ################################
    data = getAtomsInformation(filename, index)
    
    # Extract x, y, z coordinates from the data
    x2 = data[:, 2]
    y2 = data[:, 3]
    z2 = data[:, 4]

    # Set radius
    radius = 4

    for i in range(len(x2)):
        ax.scatter(x2[i], y2[i], z2[i], s=radius, c='red')

    # 关闭坐标轴网格线
    ax.grid(False)
    
    # 设置坐标轴范围
    ax.set_xlim([min(X), max(X)])
    ax.set_ylim([min(Y), max(Y)])
    ax.set_zlim([min(Z), max(Z)])
    
    # 设置按照实际数值显示比例
    ax.axis('equal')
    
    # 隐藏坐标轴
    ax.axis('off')
    
    # Show the plot
    plt.show()

others

输出起始和终止时间步

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def beginAndStoptimeStep(filename):
    target_start = "ITEM: TIMESTEP"  # 目标字符串

    first_next_line = None  # 第一次出现目标字符串的下一行内容
    last_next_line = None  # 最后一次出现目标字符串的下一行内容
    
    with open(filename, 'r') as file:
        lines = file.readlines()
        start_index = None
        end_index = None
    
        for i, line in enumerate(lines):
            if target_start in line:
                if start_index is None:
                    start_index = i + 1
                end_index = i + 1
    
        if start_index and end_index:
            first_next_line = lines[start_index].strip()
            last_next_line = lines[end_index].strip()
    
    return first_next_line, last_next_line