利用分子动力学软件LAMMPS计算固体表面对分子的吸附可以通过计算势能平均力 (Potential of Mean Force, PMF)与吸附距离的关系实现。PMF的计算一般可通过自由采样和伞形采样实现。由于伞形采样是一种常用且较为准确的方法，因此本文以NaCl表面对C2H6O吸附能的计算介绍一下伞形采样的使用。
LAMMPS in文件如下(init.data文件在文末)：
variable sigma equal 3.405 # Angstromvariable epsilon equal 0.238 # Kcal/molvariable U0 equal 10*${epsilon} # Kcal/molvariable dlt equal 0.5 # Angstromvariable x0 equal 5.0 # Angstromvariable k equal 1.5 # Kcal/mol/Angstrom^2units real # style of units (A, fs, Kcal/mol)atom_style full # molecular + chargeboundary p p p # periodic boundary conditionsread_data init.datapair_style lj/cut/coul/long 10 # cut-off 1 nmpair_coeff 1 1 0.0552 2.31pair_coeff 2 2 0.1004 4.3pair_coeff 3 3 0.0663 3.5812pair_coeff 4 4 0.0281 2.3734pair_coeff 5 5 0.0000 0.0000pair_coeff 6 6 0.2450 2.8114pair_coeff 7 7 0.1195 3.1000bond_style harmonicbond_coeff 1 349.274 1.0900bond_coeff 2 399.792 1.4300bond_coeff 3 299.844 1.5200bond_coeff 4 349.928 1.1000bond_coeff 5 442.161 0.9720angle_style harmonicangle_coeff 1 47.694 109.500angle_coeff 2 55.127 109.500angle_coeff 3 55.114 111.000angle_coeff 4 52.477 107.570angle_coeff 5 47.502 108.530angle_coeff 6 54.993 110.300dihedral_style harmonicdihedral_coeff 1 0.30115 1 3dihedral_coeff 2 1.415 1 3dihedral_coeff 3 0.0 1 1kspace_style pppm 1.0e-4pair_modify mix arithmetic tail yesneigh_modify every 1 delay 4 check yesgroup topull type 3 4 5 6 7minimize 1e-4 1e-6 100 1000reset_timestep 0#variable U atom ${U0}*atan((x+${x0})/${dlt}) &-${U0}*atan((x-${x0})/${dlt})#variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt} &-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}#fix pot all addforce 0.0 0.0 v_F energy v_Ufix mynve all nvefix mylgv all langevin 119.8 119.8 50 1530917timestep 2.0thermo 100000run 500000reset_timestep 0dump mydmp all atom 1000000 dump.lammpstrjvariable a loop 50label loopvariable zdes equal ${a}-25variable zave equal zcm(topull,z)fix mytth topull spring tether ${k} 0 0 ${zdes} 0run 200000fix myat1 all ave/time 10 10 100 v_zave v_zdes &file data-k1.5/position.${a}.datrun 1000000unfix myat1next ajump SELF loop
运行完成后需要用伞形抽样后处理软件wham进行后处理
wham下载地址: https://lammpstutorials.github.io/lammpstutorials-inputs/level3/free-energy-calculation/BiasedSampling/wham下载后，执行chmod 777 wham即可运行wham.
运行wham之前，首先使用如下python脚本处理LAMMPS运行结果：
import osk=1.5folder='data-k1.5/'f = open("metadata.dat", "w")for n in range(-50,50):datafile=folder+'position.'+str(n)+'.dat'if os.path.exists(datafile):# read the imposed position is the expected onewith open(datafile) as g:_ = g.readline()_ = g.readline()firstline = g.readline()imposed_position = firstline.split(' ')[-1][:-1]# write one file per filef.write(datafile+' '+str(imposed_position)+' '+str(k)+'\n')f.close()
运行wham
./wham -25 25 50 1e-8 119.8 0 metadata.dat PMF.dat将NaCl表面设置为0可得：
若有不足之处欢迎批评指正，
祝大家科研愉快！！
参考：
https://lammpstutorials.github.io/sphinx/build/html/tutorials/level3/free-energy-calculation.html
