lammps 怎么编译gpu版本_关于Lammps计算声子态密度（PDOS）的Matlab接口文件改进

@泛柏舟 的文章详细介绍了lammps计算声子态密度（PDOS）的几种方法，其中使用lammps输出原子速度轨迹，再采用樊哲勇老师所写matlab代码进行计算是相对可取的一种方法。

泛柏舟：通过lammps输出速度计算VODS​zhuanlan.zhihu.com

其中比较重要的部分是将lammps速度轨迹文件转换为樊哲勇matlab代码的输入文件，@泛柏舟 博文中也给出了相应的matlab代码，本文基于此版本进行了简化和改进，提升了其运算速度。

1.首先同样在Lammps中输出原子速度轨迹文件v_output.lammpstrj

只需在lammps代码中加入这样一行命令

dump 1 all custom 1 v_output.lammpstrj id type vx vy vz

2.利用matlab代码，将v_output.lammpstrj文件转换为v_out.txt文件，使其可以直接输入到樊哲勇代码中来

matlab代码如下：

function Rebuild_data(Nf,N)
% 将lammps输出的v_output.lammpstrj,转换为樊哲勇plot_pdos.m的输入文件
% Nf: 计算总时间步数
% N: 原子总数

try
    dump = fopen('v_output.lammpstrj','r'); % 读取数据文件
catch
    error('Dumpfile not found!');
end
fid = fopen('v_out.txt','w'); % 输出的后处理文件
 
for i = 1 : Nf+1 
    
    % lammps里跑了N步，实际上跑了N+1步，从0开始计算
    % 跳过前9行读取数据，注意：每次读完文件指针都不会
    % 返回初始位置，而是停留在当前位置，所以不能改变9
    
    A = textscan(dump,'%f %f %f %f %f','headerlines',9);
    data = zeros(N,5); % 储存每个时间步下的原子速度,[id type,vx,vy,yz]
    
    for j = 1 : 5
        data(:,j) = A{1,j}; % 将读取的数据（cell）存储为矩阵
    end
    
    % 重新按照升序排列的原子信息，fprintf是按列存储，所以需转置
    data = sortrows(data,1)';  
    fprintf(fid,'%d %d %f %f %fn',data); % 将原子信息存储到v_out,txt中
end

3.最后将樊哲勇matlab程序 plot_pdos.m 改写为以下命令即可

clear; close all;

% set up some parameters related to the MD simulation
N = 7480;   % number of atoms
Nc = 1000;   % number of correlation steps (a larger number gives a finer resolution)
dt = 0.001; % time interval in units of ps (its inverse is roughly the maximum frequency attainable)
omega = 1 : 0.5 : 380; % angular frequency in units of THz
nu = omega / 2 / pi; % omega = 2 * pi * nu 计算频率范围
Nf = Nc*10;

% 将lammps 输出的速度文件转换为该程序输入文件所需格式
Rebuild_data(Nf,N);

% load and rearrange the data
data = load('v_out.txt');    % get data from your output file 
Nf = size(data, 1) / N;      % number of frames
v_all = zeros(N, 3, Nf);     % all the velocity data
for n = 1 : Nf               % put the data into a good form
   v_all(:, :, n) = data((n - 1) * N + 1 : n * N, 3:5);
end