lammps运行命令行中输入变量并在MATLAB中读取log文件绘图

---

使用lammps开展模拟时，有时候我们仅想修改某一个变量并对比该变量对计算结果的影响，为了避免频繁的修改in文件，可以使用variable loop、next、jump等命令开展循环计算，也可以通过matlab循环向lammps运行命令行中输入in文件中需要修改的参数，前一种方法在之前写的文章中都已有应用，详细介绍可以参考博客《Lammps将当前体系温度作为输入，循环升温》，本文将主要介绍第二种方法，在输入参数的同时读取log文件目标参数进行绘图。

使用Al99.eam.alloy力场，计算不同晶格常数(3-5Å)下原子的平均势能，最后结果显示在4Å左右原子平均势能最小，接近Al的真是晶格常数。

---

1、说明部分

(1) in文件中参数输出

variable部分就是in文件中正常的变量设置，需要说明的是print中最后一行中%%的作用是进行标记，便于在MATLAB中找到ecoh变量。此外，最后一样print中的内容去掉%%后也是matlab代码的形式，需要在matlab中使用eval函数运行该语句来更新ecoh变量。

variable natoms equal "count(all)" 
variable teng equal "pe"
variable ecoh equal "v_teng/v_natoms"
variable L equal "lx"

print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms)= ${L};"   								
print "Cohesive energy (eV) = ${ecoh};"
print "%%  ecoh = ${ecoh};"

(2) MATLAB中循环修改in文件变量

在in文件中，晶格常数采用的变量形式(${latconst})，在matlab中设置lammps运行命令，latconst变量使用循环的方式更改。

for i = 3.0:0.10:5.0
    command_line = ['lmp_serial -var latconst ' num2str(i) ' < calc_fcc_ver2.in'];
    system(command_line)
end

(3) MATLAB中变量查找

if判断是找到%%且在该行首的行，也就是ecoh；
teval表示print的内容，即ecoh = -2.36885241762669;
eval表示执行teval中的matlab代码，也就是读取了ecoh

fid = fopen('log.lammps');
	tline = fgetl(fid);
	while ~feof(fid)
		matches = strfind(tline, '%%');
		num = length(matches);
		if num > 0 && matches == 1                                       
			teval = strrep(tline,'%%','');                              
			eval(teval)                                                 
		end
		tline = fgetl(fid);
	end
fclose(fid);

---

$附上\mathrm{i}\mathrm{n}文件及\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{l}\mathrm{a}\mathrm{b}代码$

2、in文件

# %%writefile calc_fcc_ver2.in
######################################
# LAMMPS INPUT SCRIPT
# Find minimum energy fcc configuration
# Mark Tschopp
# This requires the variable latconst to be input via the command line
# e.g., lmp_serial -var latconst 4 < calc_fcc_ver2.in

######################################
# INITIALIZATION
clear 
units metal 
dimension 3 
boundary p p p 
atom_style atomic 
atom_modify map array

######################################
# ATOM DEFINITION
lattice fcc ${latconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  
region box block 0 1 0 1 0 1 units lattice
create_box 1 box
create_atoms 1 box
replicate 1 1 1

######################################
# DEFINE INTERATOMIC POTENTIAL
pair_style eam/alloy 
pair_coeff * * Al99.eam.alloy Al
neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
######################################
# RUN 0
run 0

variable natoms equal "count(all)" 
variable teng equal "pe"
variable ecoh equal "v_teng/v_natoms"
variable L equal "lx"

print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms)= ${L};"   								
print "Cohesive energy (eV) = ${ecoh};"
print "%%  ecoh = ${ecoh};"

print "All done!" 

3、matlab代码

% MATLAB Script for running LAMMPS multiple times
clc;clear;
t_start_all = cputime;
count = 0;
for i = 3.0:0.10:5.0
    command_line = ['lmp_serial -var latconst ' num2str(i) ' < calc_fcc_ver2.in'];

    % this next line executes the command line
    system(command_line)

%   % all that is left is to mine the 'log.lammps' file for the energy
    fid = fopen('log.lammps');
        tline = fgetl(fid);
        while ~feof(fid)
           matches = strfind(tline, '%%');
           num = length(matches);
           if num > 0 && matches == 1                                       % 找到%%且在该行首
                teval = strrep(tline,'%%','');                              % 此时teval表示ecoh = -2.36885241762669;
                eval(teval)                                                 % 执行teval中的matlab代码，也就是读取了ecoh
           end
           tline = fgetl(fid);
        end
    fclose(fid);

   % store the values in a matrix
   count = count + 1;
   X(count) = i; Y(count) = ecoh;

end

plot(X,Y,'-*r'), axis square

t_start_stop = (cputime-t_start_all)/60;
fprintf('Now the used time is: %.1f mins.\n',t_start_stop);