单位

表面张力单位常用的有好几个,$\;$

• 国际制单位: Pa-m, N/m, J/m^2

• 传统单位: dyn/cm

• GROMACS单位: bar-nm

• 力场单位: kJ/mol-nm^2, kcal/mol-A^2

换算关系

1 dyn/cm = 10^-3^ Pa-m = 1 mJ/m^2

1 bar-nm = 10^-4^ Pa-m = 10^-4^ N/m = 10^-4^ J/m^2 = 0.1 mN/m = 0.1 mJ/m^2 = 0.1 dyn/cm

计算

GROMACS 中给出的 #SurfTen 为表面张力乘上表面数目, 一般slab构型有两个表面, 所以实际表面张力为

γ=#SurfTen2bar-nm=#SurfTen20mJ/m2

长程校正

对LJ势能函数, 其参数为 σ,ε, 截断距离为 Rc, 则其表面张力的长程校正为

其中

ρ(ξ) 为分子z方向的数密度分布.

这一校正项可利用matlab的二重积分函数quad2d计算, 但由于边界的不规整性, 有时不易收敛. 更好的办法是解析求解或分区域数值积分, 但实现麻烦些. 另外, 一般得到的数密度分布为离散值, 而积分需要连续值, 可以先将离散的数密度分布用样条函数进行插值. 插值方法可参考前面的一篇博文样条函数插值拟合.

测试用例

文件Rz.xvg, 单位为nm, 1/nm^3

  0                00.12      6.28082e-050.24      6.28082e-050.36      6.28082e-050.48      6.28082e-05 0.6                00.72      6.28082e-050.84      6.28082e-050.96      6.28082e-051.08      6.28082e-05 1.2      6.28082e-051.32      6.28082e-051.44      6.28082e-051.56                01.68      6.28082e-05 1.8      6.28082e-051.92      6.28082e-052.04      6.28082e-052.16      6.28082e-052.28                0 2.4      6.28082e-052.52      6.28082e-052.64      6.28082e-052.76      6.28082e-052.88      6.28082e-05   3      6.28082e-053.12                03.24      6.28082e-053.36      0.0007536983.48        0.0243696 3.6         0.3293033.72          2.270143.84           8.85253.96          20.45864.08          30.5037 4.2          33.11484.32          33.29544.44          33.15134.56          33.12534.68          32.7827 4.8          32.94094.92          33.08595.04          32.86745.16          33.06885.28          33.0544 5.4           32.9145.52           32.9565.64          33.11355.76          33.02225.88          32.9612   6          32.99316.12          33.04416.24          32.93336.36           32.9216.48          33.0009 6.6          32.91776.72            32.976.84          32.93866.96          33.15597.08           32.909 7.2          32.86647.32          33.04447.44          33.28587.56          33.34427.68          33.0472 7.8           30.7127.92          21.44868.04          9.266168.16          2.510078.28         0.383444 8.4        0.03360248.52       0.001381788.64      6.28082e-058.76      6.28082e-058.88      6.28082e-05   9      6.28082e-059.12                09.24      6.28082e-059.36      6.28082e-059.48      6.28082e-05 9.6      6.28082e-059.72      6.28082e-059.84                09.96      6.28082e-0510.08      6.28082e-0510.2      6.28082e-0510.32      6.28082e-0510.44                010.56      6.28082e-0510.68      6.28082e-0510.8      6.28082e-0510.92      6.28082e-0511.04      6.28082e-0511.16                011.28      6.28082e-0511.4      6.28082e-0511.52      6.28082e-0511.64      6.28082e-0511.76      6.28082e-0511.88      6.28082e-05

得到 Δγlrc=0.0282kcal/mol-A2=19.592mJ/m2

代码# Language: matlabfunction GammaLRC    clc;    global eps sgm Rc A B D E pp    sgm =3.15;% angstrom, TIP4P    eps =0.185207656385689;% kcal/mol =93.2kb    Rc  =2.5*sgm;    S =(sgm/Rc)^6;    A =6*pi*eps*S*(S-1);    B =3*pi*eps*Rc^2*S*(1-0.8*S);    D =3.6*pi*eps*sgm^12;    E =-3*pi*eps*sgm^6;[Z, R]=textread('Rz.xvg');    Z = Z*10;    R = R/1000;    Zmin =min(Z(:));    Zmax =max(Z(:));    pp =csape(Z,R,'second',[0,0]);%     xsp=Zmin:0.01:Zmax;%     ysp=ppval(pp,xsp);%plot(Z,R,'o',xsp,ysp,'-')    Glrc =quad2d(@PI, Zmin,Zmax, Zmin,Zmax,...'abstol',1e-6,'reltol',1e-9,'maxfunevals',8000,...'FailurePlot',true,'Singular',true)end%%function dPi =PI(z1, z2)    global Rc A B D E pp    Rz1 =ppval(pp, z1);    Rz2 =ppval(pp, z2);    ksi =abs(z1-z2);    dPi =zeros(size(ksi));    idx =(find(ksi<=Rc));dPi(idx)=(A*ksi(idx).^2+B).*Rz1(idx).*Rz2(idx);    idx =(find(ksi>Rc));dPi(idx)=( D*ksi(idx).^(-10)+E*ksi(idx).^(-4)).*Rz1(idx).*Rz2(idx);end

◆图片/表格/公式/代码完整版请参看：GROMACS表面张力单位,计算及其长程校正◆