基于 XDATCAR创建RMC_POT初始.cfg构型

1. RMC模拟与其初始构型.cfg文件

逆蒙特卡洛(RMC)模拟是根据一维衍射数据（X-ray, Neutron）反推出三维原子排列的有效手段，该算法的实现可以借助Orsolya Gereben等人开发的RMC_POT程序包。RMC_POT做RMC模拟所需的文件包括：主程序文件(.dat)，初始构型文件(.cfg)，实验数据文件(格式用户自己定义)，sfactorcube文件(控制原子移动最大位移)。在这其中.cfg初始构型文件是极其重要的。以反推金属玻璃(Metallic glasses，MG)中原子排布为例，总的来说.cfg文件的创建需要满足以下几个原则：i) 构型内原子的数密度需要符合实验测定的材料密度；ii) 原子无规堆积，缺少长程有序结构； iii) 原子两两之间的最小间隔满足定义的cutoff。如何事半功倍的构建出一组.cfg初始构型文件对进行后续独立的RMC模拟实验是至关重要的。

.cfg文件的格式(以四元MG PdNiCuP为例）如下所示：

(Version 3 format configuration file) !file created by SimpleCfg::save !
MG_Pd43Ni20Cu27P10 metallic glass@RMC_POT

0         0         0         0 moves generated, tried, accepted, potential-accepted
          0                      configurations saved

        5400 molecules of all types
          4 types of molecule
          1 is the largest number of atoms in a molecule
          0 Euler angles are provided

          F (box is cubic)
            Defining vectors are:
              20.815651   0.000000   0.000000
              0.000000   20.815651   0.000000
              0.000000   0.000000   20.815651

2322 molecules of type 1
,          1 atomic sites
              0.000000   0.000000   0.000000

1080 molecules of type 2
,          1 atomic sites
              0.000000   0.000000   0.000000

1458 molecules of type 3
,          1 atomic sites
              0.000000   0.000000   0.000000

540 molecules of type 4
,          1 atomic sites
              0.000000   0.000000   0.000000

-0.868839994050115 -0.991980000410556 -0.436619997026109
-0.808699995279919 -0.878679998204117 -0.613059997554325
-0.760199993841675 -0.952720001332732 -0.487940013432776

……

每个原子的初始位置以（-1,1）fraction坐标表示，构型满足周期性边界条件约束。

2. AIMD模拟和XDATCAR文件

AIMD模拟不依赖于势函数的限制，因而其适合处理多组元MG的模拟问题。但是AIMD建模MG存在理论及技术上的限制：a) 从玻璃弛豫现象上来讲，金属玻璃在低温状态处于绝对的亚稳定过程，弛豫时间近似无限长，理论上任何的计算机模拟/实验研究都不能遍历如此缓慢的动力学过程；b) AIMD模拟基于求解薛定谔方程，其运算速度慢，所需计算核心多，很难完整勾勒出MG实际形成过程的结构重组行为，因而基于melting-quenching创建的MG构型严格而言与实验得到的结构存在不小的差异。基于以上所述，基于AIMD输出的构型，利用衍射实验的数据做RMC精修，是捕捉实际MG玻璃结构的有效手段。

XDATCAR 文件是VASP进行AIMD计算输出的每模拟步的原子坐标集合。AIMD输出的每步构型中原子排列满足.cfg文件构建的三个重要原则，但是我们并不能直接将XDATCAR中的构型拿来进行RMC模拟，最关键的原因在于AIMD模拟构型非常小，一般<300原子数，而RMC模拟构型原子数在1000~10000之间，二者存在数量级差异。如何对AIMD的输出构型改造进行RMC模拟是个有趣的问题。

XDATCAR文件的格式(以四元MG PdNiCuP为例）如下所示：

Random MG from XJ Han
1.00000000000000
    13.8765888450000006    0.0000000000000000    0.0000000000000000
     0.0000000000000000   13.8765888450000006    0.0000000000000000
     0.0000000000000000    0.0000000000000000   13.8765888450000006
   Pd   Ni   Cu   P
    86    40    54    20
Direct
0.1967512433709550 0.0120249114084959 0.8450697621079324
0.2869488440860818 0.1819650879617962 0.5803966400872402
0.3597098725726727 0.0709252575579089 0.7680881170493263

……

每个原子的初始位置以（0,1）fraction坐标表示，构型满足周期性边界条件约束。

3. 技术路线

<1>首先取出XDATCAR中任一一步原子构型，将其用VESTA打开，另存为.cif格式文件

<2>利用materials studio打开上述.cif文件，利用build-supercell命令将其扩充(一般选择3x3x3)

<3>将上述超胞保存为.cif文件，再次利用VESTA打开，将其另存为POSCAR格式(其与上述XDATCAR格式完全一样，但是只包含一步的原子坐标)，保存时选择笛卡尔坐标

<4>利用以下的MATLAB程序段将其转变为.cfg格式

%  This code was programmed in order to transform the POSCAR data into .cfg files for RMC simulation.
% created by Ge Xuan at 2020/03/09
% Copyright interpretation belongs to Dr. Getic
format longclear;
close all;
clc;%以下循环是将模拟过程中XDATCAR构型数据依次输进程序中
fp=fopen('filepath\POSCAR_3x3x3','rt');
% c=fscanf(fp,'%c');%以字符串的形式将带题头的全部数据输入到XDATCAR矩阵中
[b1,b2,b3]=textread('filepath\POSCAR_3x3x3','%s %s %s');
XDATCAR=[b1,b2,b3];
fclose(fp);
num_atoms=input('please give the number of atoms in youe simulation box: \n');
num_steps=size(XDATCAR,1)/(num_atoms+10) % 计算出模拟步数；
if num_steps~=fix(num_steps)fprintf('you gave a wrong input')return;
endfor t=1:num_steps  %%以下的循环是按照RMC_cfg的格式输出dir = strcat('filepath\MG_PdNiCuP_', num2str(t), '.cfg');fid=fopen(dir,'wt');%在文件第一行输出固定格式字符串fprintf(fid,' %s\n','(Version 3 format configuration file) !file created by SimpleCfg::save !');%在文件第二行输出固定格式字符串和模拟步数fprintf(fid,'%s\n','MG_Pd43Ni20Cu27P10 metallic glass@RMC_POT');  %this sentence should be identified every times%第三/四行输出空白行fprintf(fid,' \n\n');%第五行输出固定的表达式fprintf(fid,'          %d         %d         %d         %d %s\n',0,0,0,0,'moves generated, tried, accepted, potential-accepted');%第六行输出固定的表达式fprintf(fid,'          %d                      %s\n \n',0,'configurations saved');%第七行输出固定的表达式fprintf(fid,'        %d %s\n',num_atoms,'molecules of all types');%第八-第十一行输出固定的表达式fprintf(fid,'          %d %s\n          %d %s\n          %d %s\n \n',4,'types of molecule',1,'is the largest number of atoms in a molecule',0,'Euler angles are provided');%第十二行到第十七行输出模拟盒子的大小fprintf(fid,'          %s\n            %s \n              %8.6f   %8.6f   %8.6f\n              %8.6f   %8.6f   %8.6f\n              %8.6f   %8.6f   %8.6f\n\n',...'F (box is cubic)','Defining vectors are:',...0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+3,1)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+3,2)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+3,3)),...0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+4,1)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+4,2)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+4,3)),...0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+5,1)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+5,2)),0.5*str2double(XDATCAR((t-1)*(num_atoms+10)+5,3)));%第十八行到第二十一行输出第一类原子的信息,这个每次都要修改一下fprintf(fid,'         %d %s\n,          %d %s \n              %7.6f   %7.6f   %7.6f\n\n',2322,'molecules of type  1',1,'atomic sites',0,0,0);%第二十二行到第二十五行输出第二类原子的信息，这个每次都要修改一下fprintf(fid,'         %d %s\n,          %d %s \n              %7.6f   %7.6f   %7.6f\n\n',1080,'molecules of type  2',1,'atomic sites',0,0,0);%第二十五行到第二十八行输出第三类原子的信息，这个每次都要修改一下fprintf(fid,'         %d %s\n,          %d %s \n              %7.6f   %7.6f   %7.6f\n\n',1458,'molecules of type  3',1,'atomic sites',0,0,0);%第二十八行到第三十一行输出第四类原子的信息，这个每次都要修改一下fprintf(fid,'         %d %s\n,          %d %s \n              %7.6f   %7.6f   %7.6f\n\n',540,'molecules of type  4',1,'atomic sites',0,0,0);%提取盒子大小x_size = str2double(XDATCAR((t-1)*(num_atoms+10)+3,1));y_size = str2double(XDATCAR((t-1)*(num_atoms+10)+4,2));z_size = str2double(XDATCAR((t-1)*(num_atoms+10)+5,3));for j=(t-1)*(num_atoms+10)+11:1:t*(num_atoms+10)fprintf(fid,' %16.15f  %16.15f  %16.15f\n',str2double(XDATCAR(j,1))/x_size*2-1,str2double(XDATCAR(j,2))/y_size*2-1,str2double(XDATCAR(j,3))/z_size*2-1);endfclose(fid);
end
toc

上述代码中的空格和换行是按照.cfg文件格式设置的，不建议随便删改；程序运行时间在4~5s之间(原子数5400个)，因而上述代码可以迅速根据XDATACR文件创造出一系列满足要求的RMC初始构型文件。

上述代码如有使用不便，请send Email to：gejade9254@163.com。谢谢