LAMMPS 中随机替换原子

发表于 更新于 分类于

在使用 LAMMPS 构建模型时，有时需要在主组分原子中随机替换为其他种类原子。比如 Fe 的 10×20×30 晶格中 15% Ni 被替换，可以通过以下方式进行替换。以下使用的脚本和文件见百度云盘，提取码 ei4i。推荐使用 LAMMPS 内建 shell+Python 命令替换。

set type/fraction 命令

set type/fraction是 LAMMPS 的内建命令，可以指定选中类型的原子一部分为另一种原子，用法如下：

1
set type type-1-ID type/fraction type-2-ID fraction seed

以上表示type-1-ID类型的原子中有fraction比例被替换成type-2-ID类型原子。比如我们实现上述 Fe 晶格中部分被替换成 Ni，LAMMPS 脚本（test.in）如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
##### 初始化（Initialization）
dimension		3
boundary		p p p
units			metal
atom_style		atomic
neighbor		0.3 bin
neigh_modify	delay 5
##### 原子参数定义（Atom definition）
lattice			bcc 2.863
region			box block 0 10 0 20 0 30
create_box		2 box
create_atoms	1 box
##### 设置（Settings）
set 			type 1 type/fraction 2 0.15 123456 # 15%的Fe被替换成Ni
mass			* 1.0 # 此处我们不关心质量
write_data		test.data

打开日志文件（test.lammps），我们发现总原子数为 12000，然而 Ni 原子数为1837，所以，set type/fraction 不能精确替换原子。

对 data 文件进行文本操作

该方式操作自由度高，甚至任意更改任意原子的属性。但对于比较大的 data 文件，通过其它软件处理缓慢。此处采用 Matlab 软件处理 data 文件。如下所示，最简单 data 文件（test1.data）格式如下：

现在我们要把Atom # atomic中部分原子的 type 替换为 2。思路如下：

  • 从 data 文件中获取总原子数
  • 设定替换比例，计算替换原子数，随机从所有原子 ID 中选取替换原子数量的 ID
  • 读取 data 文件，将上一步选取的 ID 的原子类型更改为 2

Matlab 处理脚本（test.m）如下：

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
clear; close all; clc;
C = 15/100; % 替换比例
Fid1 = fopen('test2.data', 'w+');
Fid = fopen('test1.data', 'rt');
Temp = cell(15, 1); % 前 15 行
for i = 1:15
    Temp{i} = fgets(Fid);
end
fprintf(Fid1, '%s%s%s', char(Temp(1)), char(Temp(2)), char(Temp(3)));
Temp1 = strsplit(Temp{3}, ' '); % 分割第 3行
NAtoms = str2double(Temp1{1, 1}); % 所有原子数
NReplace = round(NAtoms*C); % 替换原子数量
fprintf(Fid1, '%s\n', '2 atom types');
fprintf(Fid1, '%s%s%s%s%s%s%s%s', char(Temp(5)), char(Temp(6)), ...
    char(Temp(7)), char(Temp(8)), char(Temp(9)), char(Temp(10)), ...
    char(Temp(11)), char(Temp(12)));
fprintf(Fid1, '%s\n', '2 1');
fprintf(Fid1, '%s%s%s', char(Temp(13)), char(Temp(14)), char(Temp(15)));
rng(123456);
IDReplace = randperm(NAtoms, NReplace);
% 替换部分原子的原子类型
for i = 1:NAtoms
    Temp = fgetl(Fid);
    Temp1 = strsplit(Temp, ' ');
    if ismember(str2double(Temp1{1,1}), IDReplace) == true % 判断是否替换
        Temp1{1,2} = '2';
    end
    Temp1 = [Temp1{1,1}, ' ', Temp1{1,2}, ' ', Temp1{1,3}, ' ', ...
        Temp1{1,4}, ' ', Temp1{1,5}, ' ', Temp1{1,6}, ' ', ...
        Temp1{1,7}, ' ', Temp1{1,8}];
    fprintf(Fid1, '%s\n', Temp1);
end
Temp = fgets(Fid);
fprintf(Fid1, '%s', Temp);
Temp = fgets(Fid);
fprintf(Fid1, '%s', Temp);
Temp = fgets(Fid);
fprintf(Fid1, '%s', Temp);
% 原子速度
for i = 1:NAtoms
    Temp = fgets(Fid);
    fprintf(Fid1, '%s', Temp);
end
fclose(Fid);
fclose(Fid1);

相关的 Matlab 的命令有fopen, fgetl, fgets, fprintf, strsplit, randperm等。在替换部分原子 type 时，虽然我们可以一次性读取所有原子数据到矩阵，即Atom # atomic和Velocities之间的数据，但如果 data 文件有几百兆，则可能爆内存。建议逐行读取，虽然缓慢，但起码还可以处理文件。

调用shell命令

LAMMPS 内建shell命令，可以调用外部命令执行外部脚本，如调用 Python 程序取随机数组。由于可以把替换比例等参数直接加入 LAMMPS 脚本中，因此可以快速生成多个替换比例的 data文件，这种方式大大提高建模的效率，尤其是 data 文件超大（> 1G）时该方式是唯一准确有效的替换方式，极力推荐。思路如下：

  • 利用 LAMMPS 自带命令（lattice等）构建未替换模型
  • 设置替换比例变量（variable），计算总原子数，将替换比例和总原子数通过shell命令导出文本
  • 利用shell命令调用 Python，通过 Python 脚本将替换的原子 ID 导出文本
  • 利用shell命令将替换的原子 ID 分组（group）
  • 设定上一步分组的原子类型

LAMMPS 脚本（test3.in）如下，其中 Natoms.txt 是导出文本：

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
log				test3.lammps
##### 初始化（Initialization）
dimension		3
boundary		p p p
units			metal
atom_style		atomic
neighbor		0.3 bin
neigh_modify	delay 5
##### 原子参数定义（Atom definition）
lattice			bcc 2.863
region			box block 0 10 0 20 0 30
create_box		2 box
create_atoms	1 box
##### 设置（Settings）
variable		Natoms equal atoms
variable		N equal ${Natoms}
variable		C equal 0.15 # 替换比例
shell			touch Natoms.txt # 计算总原子数
shell			echo ${C} > Natoms.txt # 输出替换比例到文本
shell			echo ${N} >> Natoms.txt # 输出总原子数到文本
shell			python3 test.py # 调用 Python np.random.choice 随机选择原子 ID，并输出文本
shell			"sed -i 's/^/group Ni id &/g' ID.txt" # 在文本前添加 group Ni id 命令
include			ID.txt
set				group Ni type 2
mass			* 1.0 # 此处我们不关心质量
write_data		test3.data

Python 脚本（test.py）如下，该脚本根据 LAMMPS 中导出的含原子数量和替换比例的文件 ID.txt，随机选取替换原子 ID，再将原子 ID 导出文本。

1
2
3
4
5
6
7
# -*- coding: utf-8 -*-
import numpy as np
A = np.loadtxt('Natoms.txt')
np.random.seed(123456)
NReplace = int(round(A[1]*A[0])) # 替换原子数量
IDReplace = np.random.choice(int(A[1]), size=NReplace, replace=False)+1 # 替换原子 ID
np.savetxt('ID.txt', IDReplace, fmt='%u', newline=' ')