[游戏开发]巨正则蒙特卡洛—甲烷和二氧化碳的竞争吸附GCMC

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

一、What is 巨正则蒙特卡洛？

$\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}$ 即是 $\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{c}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{M}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{C}\mathrm{a}\mathrm{r}\mathrm{l}\mathrm{o}$，巨正则蒙特卡洛。$“热浴”$ 这个翻译我觉得很形象，手册中用的是$imaginaryreservoir$。比如 $\mathrm{N}\mathrm{V}\mathrm{T}$系综，当采用 $\mathrm{N}\mathrm{V}\mathrm{T}$ 系综控温时，可以理解为把模拟体系放在一个一定分子数、恒定温度的大热浴中，模拟体系的体积是恒定的、并且不能和大热浴交换分子。显然，根据热力学定律，达到平衡态的时候，模拟体系和热浴的温度相同，这就达到了控温的目的。

那么在GCMC中，模拟体系和一个具有一定化学式 $\mu$ 和一定温度的大热浴接触（这里的大热浴是一个孤立系综）。同时允许模拟体系和大热浴交换粒子数。那么对于这样的开放系统，即粒子数 $N$是可以变化的，因此化学势 $\mu$被引入了热力学的基本方程：

$dE=TdS?PdV+\mu dN$
这里的 $\mu$就是化学式 $\mathrm{C}\mathrm{h}\mathrm{e}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$.

那么整个孤立体系（大热浴A和其中的模拟体系B）显然达到热力学平衡的时候，必须满足以下条件：

1). 温度相同：$T_A=T_B$；
2). 压力相同：$P_A=P_B$；
3). 化学势相同：${\mu }_A={\mu }_B$；

$显然化学势是对粒子扩散趋势的一种量度。$

二、$\mathrm{L}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{p}\mathrm{s}$如何实现GCMC—$\mathrm{F}\mathrm{i}\mathrm{x}\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}$

1. 核心命令详解

Lammps手册中，对GCMC命令的调用是这样的：

$这里的核心命令：$

$\mathrm{N}：每多少步来调用一次，\mathrm{F}\mathrm{i}\mathrm{x}?\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}；$
$\mathrm{X}：每\mathrm{N}步进行约\mathrm{X}次的\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}交换；$
$\mathrm{M}：每\mathrm{N}步进行尝试约\mathrm{M}个\mathrm{M}\mathrm{C}移动；$
$\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}：插入原子的类型；$
$\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{d}：随机数种子；$
$\mathrm{T}：上文提到的大热浴的温度；$
$\mu ：化学势；$
$\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}：\mathrm{M}\mathrm{C}移动的最大距离；$
$\mathrm{m}\mathrm{c}\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{s}：\mathrm{M}\mathrm{C}移动中的原子比例；$
$\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{d}:\mathrm{F}\mathrm{i}\mathrm{x}?\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{d}/\mathrm{s}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{l}命令的\mathrm{I}\mathrm{D}；$
$\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{k}\mathrm{e}:\mathrm{f}\mathrm{i}\mathrm{x}?\mathrm{I}\mathrm{D},\mathrm{f}\mathrm{i}\mathrm{x}?\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{k}\mathrm{e}的命令；$
$\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}:\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}?\mathrm{I}\mathrm{D}进行\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}的区域；$
$\mathrm{f}\mathrm{u}\mathrm{g}\mathrm{a}\mathrm{c}\mathrm{i}\mathrm{t}{\mathrm{y}}_{\mathrm{c}}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}:逸度系数；$
$\mathrm{f}\mathrm{u}\mathrm{l}{\mathrm{l}}_{\mathrm{e}}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{y}:进行\mathrm{G}\mathrm{C}\mathrm{M}\mathrm{C}交换时计算体系总能量；$
$\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}:插入原子的电荷量；$
$\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}:插入原子到哪个\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}；$
$\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}:\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}：插入原子的类型；$
$\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\mathrm{I}\mathrm{D}插入到哪个\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}；$
$\mathrm{t}\mathrm{f}\mathrm{a}\mathrm{c}\_\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{t}:按比例提高/降低插入原子的温度；$
$\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}\_\mathrm{c}\mathrm{u}\mathrm{t}\mathrm{o}\mathrm{f}\mathrm{f}:删除重叠原子的最大距离；$

关于分子模板请参考《LAMMPS 中 molecule command 的分子模板》

2. 需要注意自由度

由于体系中，允许粒子数的改变，这里需要对温度进行修正，手册中也给出了解释。

3. MD+MC的设置

当体系中的温度采用$FixNVT$ 进行控温时，就是一个GCMC+MD的模拟。当然NVT的温度必须设置为和GCMC中目标温度相同。同样的自由度的问题也需要修正。

4. $\mathrm{t}\mathrm{f}\mathrm{a}\mathrm{c}\_\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{t}$ 关键词的理解

这里由于插入的原子或者分子，根据系统目标来指定，并且分子是$\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{i}\mathrm{d}$，这会导致插入原子的“温度”较低。因此，需要采用一个因子：$\mathrm{t}\mathrm{f}\mathrm{a}\mathrm{c}\_\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{t}$，来尽快平衡。Lammps 官方给出的案例中采用了这样的方法：

# gcmc
variable        tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix             mygcmc co2 gcmc 100 100 0 0 54341 ${temp} ${mu} ${disp} mol &
                co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt

# atom counts

5. 化学势 $\mu$ 的设置

根据上面的讨论我们可以看到，插入和删除原子会导致体系的化学势 $\mu$ 和压力 $p$的改变。这里手册中对于化学势的计算是：

我们通过化简得到：

$e^{\frac{\mu }{k_{b}T}}=\frac{?P{\Lambda }^3}{k_{b}T}$
显然，化学势和压力存在以上的对应关系。那么在设置的时候，可以通过直接指定压力，无需设置化学势。这里在论坛中我也找到了对应的讨论可以参考。

$注意！！！$

指定等同於指定化學勢只對單組分系統管用。單組分系統的miu 是T和p的函數。而多組分的時候miu = f (T, p, nA, nB…)

二、代码实现—GCMC甲烷和二氧化碳的竞争吸附

1. 循环不同压力的设置

#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
# loop make all_file
echo   screen
print        " ------------------------------------------------"
print        " -----------------NOW making files---------------"
print        " ------------------------------------------------"
label  out_makefile
    variable     number_out loop 1             
    variable     file_name_out index 1  
    shell        mkdir ${file_name_out}   
    
    shell        cp CH4.txt ${file_name_out}
    shell        cp CO2.txt ${file_name_out}   
    shell        cp gra_layer.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 1  

variable     loop_file_out_name index 1 
shell        cd ${loop_file_out_name}
log             ${loop_file_out_name}.log

2. 物理模型构建

###################
# initial setting #
###################
#echo              screen
units             metal
dimension         3
boundary          p p p
atom_style        full
bond_style        harmonic
angle_style       harmonic


#----------------------------------------------------------#
variable          TS        equal 0.0001         ### ${TS} timestep
variable          Tdamp     equal 100*${TS}      ### ${Tdamp} Tdamp
variable          T         equal 300            ### ${T} initial temperature
variable          Height    equal 100             ### ${Height} height between two plate

variable          w_down    equal 18             ###
variable          w_up      equal 16+${Height}   ###
#----------------------------------------------------------#
region            box   block 0 100 0 100 0 100 units box

create_box        5 box                            &
                  bond/types 2                     &
                  angle/types 2                    &
                  extra/bond/per/atom 10           &
                  extra/angle/per/atom 10          &
                  extra/special/per/atom 10
                     
#----------------------------------------------------------#
read_data        gra_layer.data         &
                    add append nocoeff  &
                    group graphene_down
                    
read_data         gra_layer.data        &
                    add append nocoeff  &
                    shift 0 0 ${w_up}   &
                    group graphene_up
#----------------------------------------------------------#

# for CO2 
bond_coeff        1       0       1.16   
angle_coeff       1       0       180 


# for CH4 
bond_coeff        2  14.35   1.09
angle_coeff       2  1.517   109.5

#----------------------------------------------------------#
# setting mass
mass              1 12   # C

mass              2 16   # water-O
mass              3 1    # water-H

mass              4 12   # CH4-C
mass              5 16   # CH4-O

# make gas region
#----------------------------------------------------------#
region           gas_block block          &
                     20           80   &
                     20           80   &
                 30 70 units box 
                 
molecule        CO2  CO2.txt toff 1
create_atoms    0  random 1 12345  gas_block  mol  CO2  12345
#----------------------------------------------------------#

# make Ch4
molecule        CH4  CH4.txt toff 3 boff 1 aoff 1
#create_atoms    0  random 1 12345  gas_block  mol  CH4  55534
#----------------------------------------------------------#
###################
#  group setting  #
###################
#----------------------------------------------------------#
group             gra       type 1
group             CO_C      type 2
group             CO_O      type 3

group             CH_C      type 4
group             CH_H      type 5

group             graphene  type 1  
group             CO2       type 2 3
group             CH4       type 4 5

group             all_gas   type 2 3 4 5 


set              type     1 charge    0.0
set              type     2 charge    0.7
set              type     3 charge   -0.35

set              type     4 charge   -0.24
set              type     5 charge    0.06
#----------------------------------------------------------#
write_data 11.data
#----------------------------------------------------------#
#      graphene-1; CO_C-2; CO_O-3; CH_C-4; CH_H-5;
#----------------------------------------------------------#

3. 力场设置

#----------------------------------------------------------#
#      graphene-1; CO_C-2; CO_O-3; CH_C-4; CH_H-5;
#----------------------------------------------------------#

pair_style        lj/cut/coul/long 12
pair_modify       mix arithmetic tail yes

pair_coeff        1 1  0.003034  3.55  # gra-gra

# rigid CO2 TraPPE model 

pair_coeff        2 2  0.053649   2.8
pair_coeff        3 3  0.156973   3.05 

pair_coeff        4 4  0.002864  3.5   # CH_C CH_C
pair_coeff        5 5  0.001300  2.5   # CH_C CH_C


neigh_modify      exclude type 1 1

#--------------------------------------------------------------#

kspace_style      pppm 1.0e-5 1.0e-6

velocity          all_gas create ${T} 12345 loop local #

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

neighbor        2 bin
neigh_modify    every 1 delay 10 check yes

4. MD+GCMC模拟

# rigid constraints with thermostat 
fix              CO2_nvt      CO2 rigid/nvt/small molecule temp ${T} ${T} ${Tdamp} mol CO2
fix_modify       CO2_nvt      dynamic/dof yes

# rigid constraints with thermostat 
#fix              CH4_nvt      CH4 rigid/nvt/small molecule temp ${T} ${T} ${Tdamp} mol CH4
#fix_modify       CH4_nvt      dynamic/dof yes


fix              mygcmc_CO2 CO2 gcmc 100 100 0 0 54341 ${T} 0 1 mol &
                 CO2 pressure 10 region box rigid CO2_nvt

fix              mygcmc_CH4 CH4 gcmc 100 100 0 0 12321 ${T} 0 1 mol &
                 CH4 pressure 10 region box 


compute_modify   thermo_temp dynamic/dof yes

三、全部代码

###################################################################
#             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 1             
    variable     file_name_out index 1  
    shell        mkdir ${file_name_out}   
    
    shell        cp CH4.txt ${file_name_out}
    shell        cp CO2.txt ${file_name_out}   
    shell        cp gra_layer.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 1  

variable     loop_file_out_name index 1 
shell        cd ${loop_file_out_name}
log             ${loop_file_out_name}.log


###################
# initial setting #
###################
#echo              screen
units             metal
dimension         3
boundary          p p p
atom_style        full
bond_style        harmonic
angle_style       harmonic


#----------------------------------------------------------#
variable          TS        equal 0.0001         ### ${TS} timestep
variable          Tdamp     equal 100*${TS}      ### ${Tdamp} Tdamp
variable          T         equal 300            ### ${T} initial temperature
variable          Height    equal 100             ### ${Height} height between two plate

variable          w_down    equal 18             ###
variable          w_up      equal 16+${Height}   ###
#----------------------------------------------------------#
region            box   block 0 100 0 100 0 100 units box

create_box        5 box                            &
                  bond/types 2                     &
                  angle/types 2                    &
                  extra/bond/per/atom 10           &
                  extra/angle/per/atom 10          &
                  extra/special/per/atom 10
                     
#----------------------------------------------------------#
read_data        gra_layer.data         &
                    add append nocoeff  &
                    group graphene_down
                    
read_data         gra_layer.data        &
                    add append nocoeff  &
                    shift 0 0 ${w_up}   &
                    group graphene_up
#----------------------------------------------------------#

# for CO2 
bond_coeff        1       0       1.16   
angle_coeff       1       0       180 


# for CH4 
bond_coeff        2  14.35   1.09
angle_coeff       2  1.517   109.5

#----------------------------------------------------------#
# setting mass
mass              1 12   # C

mass              2 16   # water-O
mass              3 1    # water-H

mass              4 12   # CH4-C
mass              5 16   # CH4-O

# make gas region
#----------------------------------------------------------#
region           gas_block block          &
                     20           80   &
                     20           80   &
                 30 70 units box 
                 
molecule        CO2  CO2.txt toff 1
create_atoms    0  random 1 12345  gas_block  mol  CO2  12345
#----------------------------------------------------------#

# make Ch4
molecule        CH4  CH4.txt toff 3 boff 1 aoff 1
#create_atoms    0  random 1 12345  gas_block  mol  CH4  55534
#----------------------------------------------------------#
###################
#  group setting  #
###################
#----------------------------------------------------------#
group             gra       type 1
group             CO_C      type 2
group             CO_O      type 3

group             CH_C      type 4
group             CH_H      type 5

group             graphene  type 1  
group             CO2       type 2 3
group             CH4       type 4 5

group             all_gas   type 2 3 4 5 


set              type     1 charge    0.0
set              type     2 charge    0.7
set              type     3 charge   -0.35

set              type     4 charge   -0.24
set              type     5 charge    0.06
#----------------------------------------------------------#
write_data 11.data
#----------------------------------------------------------#
#      graphene-1; CO_C-2; CO_O-3; CH_C-4; CH_H-5;
#----------------------------------------------------------#

pair_style        lj/cut/coul/long 12
pair_modify       mix arithmetic tail yes

pair_coeff        1 1  0.003034  3.55  # gra-gra

# rigid CO2 TraPPE model 

pair_coeff        2 2  0.053649   2.8
pair_coeff        3 3  0.156973   3.05 

pair_coeff        4 4  0.002864  3.5   # CH_C CH_C
pair_coeff        5 5  0.001300  2.5   # CH_C CH_C


neigh_modify      exclude type 1 1

#--------------------------------------------------------------#

kspace_style      pppm 1.0e-5 1.0e-6

velocity          all_gas create ${T} 12345 loop local #

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

neighbor        2 bin
neigh_modify    every 1 delay 10 check yes

#--------------------------------------------------------------# 


# rigid constraints with thermostat 
fix              CO2_nvt      CO2 rigid/nvt/small molecule temp ${T} ${T} ${Tdamp} mol CO2
fix_modify       CO2_nvt      dynamic/dof yes

# rigid constraints with thermostat 
#fix              CH4_nvt      CH4 rigid/nvt/small molecule temp ${T} ${T} ${Tdamp} mol CH4
#fix_modify       CH4_nvt      dynamic/dof yes


fix              mygcmc_CO2 CO2 gcmc 100 100 0 0 54341 ${T} 0 1 mol &
                 CO2 pressure 10 region box rigid CO2_nvt

fix              mygcmc_CH4 CH4 gcmc 100 100 0 0 12321 ${T} 0 1 mol &
                 CH4 pressure 10 region box 


compute_modify   thermo_temp dynamic/dof yes

#--------------------------------------------------------------#

thermo            1000   
thermo_style      custom step temp press 
variable          1000ps equal 100/${TS}

dump              123  all custom 500 dump.123   &
                  id type x y z vx vy vz fx fy fz
#--------------------------------------------------------------#                  
###################
#    run 1        #
###################     

run               ${1000ps}

#run              1
#--------------------------------------------------------------#

write_data        final.data


shell            cd ..
clear

next             loop_file_out_name
jump             SELF calculation_in_outfile

四、$\mathrm{L}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{p}\mathrm{s}$ 给出的案例

# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
# Rigid CO2 TraPPE model
# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
# mixtures containing alkanes, carbon dioxide and
# nitrogen AIChE J., 47,1676-1682 (2001)].

# variables available on command line

variable        mu index -8.1
variable	disp index 0.5
variable        temp index 338.0
variable        lbox index 10.0
variable        spacing index 5.0

# global model settings

units           real
atom_style      full
boundary        p p p
pair_style      lj/cut/coul/long  14 
pair_modify     mix arithmetic tail yes
kspace_style    ewald 0.0001
bond_style      harmonic
angle_style     harmonic

# box, start molecules on simple cubic lattice

lattice 	sc ${spacing}
region          box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
create_box      2 box                       &
                bond/types 1                &
                angle/types 1               &
                extra/bond/per/atom 2       &
                extra/angle/per/atom 1      &
                extra/special/per/atom 2
molecule        co2mol CO2.txt
create_atoms   	0 box mol co2mol 464563 units box
                        
# rigid CO2 TraPPE model 

pair_coeff      1   1  0.053649   2.8
pair_coeff      2   2  0.156973   3.05 
bond_coeff      1       0       1.16   
angle_coeff     1       0       180 

# masses

mass 1 12.0107 
mass 2 15.9994 

# MD settings

group           co2 type 1 2
neighbor        2.0 bin
neigh_modify    every 1 delay 10 check yes
velocity       	all create ${temp} 54654
timestep        1.0

# rigid constraints with thermostat 

fix             myrigidnvt co2 rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol

# dynamically update  fix rigid/nvt/small temperature ndof
fix_modify	myrigidnvt dynamic/dof yes

# gcmc

variable        tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
fix             mygcmc co2 gcmc 100 100 0 0 54341 ${temp} ${mu} ${disp} mol &
                co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt

# atom counts

variable 	carbon atom "type==1"
variable        oxygen atom "type==2"
group 		carbon dynamic co2 var carbon
group 	        oxygen dynamic co2 var oxygen
variable        nC equal count(carbon)
variable        nO equal count(oxygen)

# output

variable	tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable	iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable	dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable	racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)

# dynamically update default temperature ndof
compute_modify  thermo_temp dynamic/dof yes

thermo_style    custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nC v_nO
thermo          1000

# run

run             20000