Lammps压力控制之活塞控压—一个代码循环限域空间内的气体分子数

Mr. Material

已于 2023-03-05 11:01:51 修改

阅读量4.1k

点赞数 7

分类专栏： Lammps实例教程—In文件详解文章标签： matlab Lammps 活塞控压

于 2022-03-04 22:32:54 首次发布

本文链接：https://blog.csdn.net/qq_43689832/article/details/123286736

版权

Lammps实例教程—In文件详解专栏收录该内容

33 篇文章 1034 订阅 ¥19.90 ¥99.00

订阅专栏

超级会员免费看

关注 $\color{Violet} \rm Mr.material\ ,$ $\color{red}{更}$ $\color{blue}{多}$ $\color{orange}{精}$ $\color{green}{彩}$ ！

主要专栏内容包括：
$\textbf{ \underline{\dag《LAMMPS小技巧》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关安装教程、原理以及模拟小技巧（难度： $\bigstar$ ）
$\textbf{ \underline{\dag\dag《LAMMPS实例教程—In文件详解》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关物理过程模拟。（包含：热导率计算、定压比热容计算，难度： $\bigstar$ $\bigstar$ $\bigstar$ ）
$\textbf{ \underline{\dag\dag\dag《Lammps编程技巧及后处理程序技巧》：}}$ 主要介绍针对分子模拟的动力学过程（轨迹文件）进行后相关的处理分析（需要一定编程能力。难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag《分子动力学后处理集成函数—Matlab》：}}$ 主要介绍针对后处理过程中指定函数，进行包装，方便使用者直接调用（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag《SCI论文绘图—Python绘图常用模板及技巧》：}}$ 主要介绍针对处理后的数据可视化，并提供对应的绘图模板（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag\dag《分子模拟—Ovito渲染案例教程》：}}$ 主要采用 $\rm Ovito$ 软件，对 $\rm Lammps$ 生成的轨迹文件进行渲染（难度： $\bigstar$ $\bigstar$ ）。

$\color{red}{\textbf{ \underline{专栏说明（订阅后可浏览对应专栏全部博文）：}}}$
$\color{red} 注意：$ 如需只订阅某个单独博文，请联系博主邮箱咨询。 $\rm lammps\_materials@163.com$

$\spadesuit$ $\dag$ 开源后处理集成程序：请关注专栏《LAMMPS后处理——MATLAB子函数合集整理》
$\spadesuit$ $\dag$ $\dag$ 需要付费定制后处理程序请邮件联系： $\rm lammps\_materials@163.com$

专栏浏览《LAMMPS实例教程—In文件详解》

请添加图片描述

一、LAMMPS是如何控压的

如果这里我们只考虑理想气体对压力的贡献，即： $\color{red}PV=nk_{b}T$ ，当体系内原子数不变时，不同压力是可以通过压缩盒子大小实现的。
在这里插入图片描述
$\color{red}那么这样的话，如果对于限域空间的压力控制，如何实现呢？$
$\color{red}也就是我们只想压缩气体部分，不想压缩固体部分。$

二、活塞控压—压力如何理解？

这里我们可以理解为，采用一个刚体平板（平板内原子没有相对运动）去压气体。那么这个压力，可以表示为： $\color{red}p=F/S$ 。相同的，每一个原子受力等于 $\color{blue}F_{atom}=F/all\ atom$ 。
在这里插入图片描述

三、程序如何实现刚性板的压力控制

显然，根据上面的分析，这个物理过程可以分为下面的程序：
$\color{red} 1) 刚性二氧化硅板原子在x, y方向，不受力；$
$\color{red} 2) 刚性二氧化硅板原子与气体的相互作用要计算，每一次fix我都把刚性板的力清零；$
$\color{red} 3) 刚性二氧化硅板原子受力为根据目标压强计算好的恒力；$

$\color{blue} 请注意：这里气体的压力达到平衡的时候，刚性板将会维持在一定位置；$

fix               platefix1               silica_hot  setforce 0.0 0.0 NULL
fix               platefix2               silica_hot  aveforce 0.0 0.0 0.0
fix               platefix3               silica_hot  addforce 0.0 0.0 -v_force_values
fix               NVE_silica_hot          silica_hot  nve

四、根据理想气体方程对活塞控压进行验证

这里我计算出 $100 n s$ ，每一帧的压力，可以明显看到，达到平衡后，压力和参考值符合很好。

请添加图片描述

五、全部代码实现

$\color{red}注意这里指提供编程思路，需要根据实际体系进行修改！$

###################################################################
#             Author Email: lammps_materials@163.com              #
###################################################################
#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
# loop make all_file
echo   screen
print        " ------------------------------------------------"
print        " -----------------NOW making files---------------"
print        " ------------------------------------------------"
label  out_makefile
	variable     number_out loop 2             
	variable     file_name_out index 5000 6000 
	shell        mkdir ${file_name_out}   
	
	shell        cp newsio.tersoff ${file_name_out}
	shell        cp 123.data ${file_name_out}	
	
	next         number_out
	next         file_name_out
	
jump         SELF  out_makefile	
clear
print        "--------------------------------------------"
print        "------------------ Finish! -----------------"
print        "--------------------------------------------"
#-----------------------------------------------------------------#
label        calculation_in_outfile
variable     loop_number_out loop 2  

variable     loop_file_out_name index 5000 6000
shell        cd ${loop_file_out_name}
log             ${loop_file_out_name}.log


###################
# initial setting #
###################
echo              screen
units             metal
dimension         3
boundary          p p s
atom_style        full
#----------------------------------------------------------#
variable          TS        equal 0.005         ### ${TS} timestep
variable          Tdamp     equal 100*${TS}      ### ${Tdamp} Tdamp
variable          T         equal 300            ### ${T} initial temperature
variable          Height    equal 30             ### ${Height} height between two plate

variable          w_down    equal 25             ###
variable          w_up      equal 23+${Height}   ###
#----------------------------------------------------------#
region            box   block 0 1 0 1 0 1 units box

create_box        3 box                           
#----------------------------------------------------------#
read_data         123.data              &
                    add append nocoeff  &
					group silica_cold
					
read_data         123.data              &
                    add append nocoeff  &
					shift 0 0 ${w_up}   &
					group silica_hot
#----------------------------------------------------------#
variable         x_length equal lx
variable         y_length equal ly
variable         seed equal  12345

region           water_block block        &
                     0      ${x_length}   & 
				     0      ${y_length}   &
				  ${w_down} ${w_up} units box 

lattice          sc  4.0
create_atoms     3   random ${loop_file_out_name} ${seed} water_block

#----------------------------------------------------------#
###################
#  group setting  #
###################
#----------------------------------------------------------#
# setting mass
mass              1 28.0
mass              2 28.0 # 16
mass              3 28.0
#----------------------------------------------------------#
group             si_si   type 1
group             si_o    type 2
group             N2      type 3

group             silica  type 1 2 
group             N2      type 3 

#----------------------------------------------------------#
#minimize          1.0e-5 1.0e-7 1000 10000
#----------------------------------------------------------#

#--------------------------------------------------------------#
pair_style        hybrid tersoff lj/cut 10
pair_coeff        * * tersoff newsio.tersoff Si O NULL  

### In N2  
#--------------------------------------------------------------#

pair_coeff        3 3  lj/cut  0.0082037  3.75    ### N2 and N2

#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
#       in silica surface 1 si 2 o ---> 3 o 4 h
#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
variable         scale     equal     1
#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
variable         ev_23     equal   0.012751068 
variable          A_23     equal   3.225
#---------------------------------------------------------------------------#

#variable          si_N13   equal  ${scale}*${ev_13}
variable         sio_N23   equal  ${scale}*${ev_23}
#---------------------------------------------------------------------------#
pair_coeff        1 3  lj/cut      0 0                      #  sio-N
pair_coeff        2 3  lj/cut      ${sio_N23} ${A_23}       #  sio-N

#--------------------------------------------------------------#
#--------------------------------------------------------------#
timestep          ${TS}


velocity          N2 create ${T} 12345 loop local #
minimize          1.0e-5 1.0e-7 1000 10000
neighbor          0.5 bin
neigh_modify      every 5 delay 0 check no

neigh_modify     exclude type 1 2
neigh_modify     exclude type 2 2
neigh_modify     exclude type 1 1

write_data       1.data
#--------------------------------------------------------------# 
variable          each_atom               equal    1.49*1e-5

compute           myTemp                  N2 temp
compute           center_silica           silica_hot com
fix               NVT_N2                  N2  nvt temp ${T} ${T} ${Tdamp} 
fix_modify        NVT_N2 temp myTemp 
#fix              NVT_silica              silica_cold nvt temp ${T} ${T} ${Tdamp}
#fix              silica_fix              silica spring/self ${loop_file_out_name}
#fix              silica_recenter_cold    silica_cold recenter INIT INIT INIT units box
#--------------------------------------------------------------#

fix               platefix1               silica_hot  setforce 0.0 0.0 NULL
fix               platefix2               silica_hot  aveforce 0.0 0.0 0.0
fix               platefix3               silica_hot  addforce 0.0 0.0 -v_force_values
fix               NVE_silica_hot          silica_hot  nve
thermo            10000     
thermo_style      custom step c_myTemp press vol
variable          1000ps equal 100000/${TS}
#100  100  10000
fix               55 silica_hot ave/time  100  1000  100000  c_center_silica[*]  file com.dat
#10000
dump              123  all custom 50000 dump.123   &
                  id type x y z vx vy vz fx fy fz
#--------------------------------------------------------------#				  
###################
#    run 3        #
###################				  
run               ${1000ps}
#run              10
#--------------------------------------------------------------#

write_data        final.data
write_restart     heat_conductivity.file

shell            cd ..
clear

next             loop_file_out_name
jump             SELF calculation_in_outfile

六、Lammps中的相似案例

$\color{red}注意这里指提供编程思路，需要根据实际体系进行修改！$

请添加图片描述

# 2-d LJ flow simulation

dimension	2
boundary	p s p

atom_style	atomic
neighbor	0.3 bin
neigh_modify	delay 5

# create geometry

lattice		hex 0.7
region		box block 0 20 0 10 -0.25 0.25
create_box	3 box
create_atoms	1 box

mass		1 1.0
mass		2 1.0
mass		3 1.0

# LJ potentials

pair_style	lj/cut 1.12246
pair_coeff	* * 1.0 1.0 1.12246

# define groups

region	     1 block INF INF INF 1.25 INF INF
group	     lower region 1
region	     2 block INF INF 8.75 INF INF INF
group	     upper region 2
group	     boundary union lower upper
group	     flow subtract all boundary

set	     group lower type 2
set	     group upper type 3

# initial velocities

compute	     mobile flow temp
velocity     flow create 1.0 482748 temp mobile
fix	     1 all nve
fix	     2 flow temp/rescale 200 1.0 1.0 0.02 1.0
fix_modify   2 temp mobile

# Couette flow

#velocity     lower set 0.0 0.0 0.0
#velocity     upper set 3.0 0.0 0.0
#fix	     3 boundary setforce 0.0 0.0 0.0
#fix	     4 all enforce2d

# Poiseuille flow

velocity     boundary set 0.0 0.0 0.0
fix	     3 lower setforce 0.0 0.0 0.0
fix	     4 upper setforce 0.0 NULL 0.0
fix	     5 upper aveforce 0.0 -1.0 0.0
fix	     6 flow addforce 0.5 0.0 0.0
fix	     7 all enforce2d

# Run

timestep	0.003
thermo		500
thermo_modify	temp mobile

dump		1 all atom 500 dump.flow

#dump		2 all image 100 image.*.jpg type type &
#		zoom 1.6 adiam 1.5
#dump_modify	2 pad 5

#dump		3 all movie 100 movie.mpg type type &
#		zoom 1.6 adiam 1.5
#dump_modify	3 pad 5

run		10000