2014-09-24 13:00:48
单位
表面张力单位常用的有好几个,
- 国际制单位: 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
| Unit | Pa-m N/m J/m2 | mN/m mJ/m2 dyn/cm | bar-nm | kcal/mol-A2 | kJ/mol-nm2 |
|---|---|---|---|---|---|
| Pa-m N/m J/m2 | 1 | 103 | 104 | 1.43932643164436 | 602.214178999998 |
| mN/m mJ/m2 dyn/cm | 10-3 | 1 | 10 | 1.43932643164436E3 | 602.214178999998E3 |
| bar-nm | 10-4 | 0.1 | 1 | 1.43932643164436E4 | 602.214178999998E4 |
| kcal/mol-A2 | 0.694769426875285 | 0.694769426875285E3 | 0.694769426875285E4 | 1 | 418.4 |
| kJ/mol-nm2 | 1.66053878316273E-3 | 1.66053878316273 | 16.6053878316273 | 2.39005736137667E-3 | 1 |
计算
GROMACS 中给出的 #SurfTen 为表面张力乘上表面数目, 一般slab构型有两个表面, 所以实际表面张力为
$\g = {\text{#SurfTen} \over 2} \; \text{bar-nm} ={\text{#SurfTen} \over 20} \; \text{mJ/m}^2$
长程校正
对LJ势能函数, 其参数为 $\s, \ve$, 截断距离为 $R_c$, 则其表面张力的长程校正为
\[\begin{align} \D \g^{\text{lrc} } &= \iint \r(z) \r(z') \P(\x) dz dz' \\ \P(\x) &= \p_{zz}(\x)-\p_{xx}(\x) \\ \x &=|z-z'| \end{align}\]其中
\[\P(\x) = \begin{cases} 6\p\ve \x^2[ ({\s \over R_c})^{12} - ({\s \over R_c})^6 ] + 3\p\ve R_c^2[ ({\s \over R_c})^6 - {4 \over 5} ({\s \over R_c})^6 ] \; & \x \le R_c \\ 3\p\ve \x^2[ {6 \over 5} ({\s \over \x})^{12} - ({\s \over \x})^6] \; & \x > R_c \end{cases}\]$\r(\xi)$ 为分子z方向的数密度分布.
这一校正项可利用matlab的二重积分函数quad2d计算, 但由于边界的不规整性, 有时不易收敛. 更好的办法是解析求解或分区域数值积分, 但实现麻烦些.
另外, 一般得到的数密度分布为离散值, 而积分需要连续值, 可以先将离散的数密度分布用样条函数进行插值.
插值方法可参考前面的一篇博文样条函数插值拟合.
测试用例
文件Rz.xvg, 单位为nm, 1/nm^3
0 0
0.12 6.28082e-05
0.24 6.28082e-05
0.36 6.28082e-05
0.48 6.28082e-05
0.6 0
0.72 6.28082e-05
0.84 6.28082e-05
0.96 6.28082e-05
1.08 6.28082e-05
1.2 6.28082e-05
1.32 6.28082e-05
1.44 6.28082e-05
1.56 0
1.68 6.28082e-05
1.8 6.28082e-05
1.92 6.28082e-05
2.04 6.28082e-05
2.16 6.28082e-05
2.28 0
2.4 6.28082e-05
2.52 6.28082e-05
2.64 6.28082e-05
2.76 6.28082e-05
2.88 6.28082e-05
3 6.28082e-05
3.12 0
3.24 6.28082e-05
3.36 0.000753698
3.48 0.0243696
3.6 0.329303
3.72 2.27014
3.84 8.8525
3.96 20.4586
4.08 30.5037
4.2 33.1148
4.32 33.2954
4.44 33.1513
4.56 33.1253
4.68 32.7827
4.8 32.9409
4.92 33.0859
5.04 32.8674
5.16 33.0688
5.28 33.0544
5.4 32.914
5.52 32.956
5.64 33.1135
5.76 33.0222
5.88 32.9612
6 32.9931
6.12 33.0441
6.24 32.9333
6.36 32.921
6.48 33.0009
6.6 32.9177
6.72 32.97
6.84 32.9386
6.96 33.1559
7.08 32.909
7.2 32.8664
7.32 33.0444
7.44 33.2858
7.56 33.3442
7.68 33.0472
7.8 30.712
7.92 21.4486
8.04 9.26616
8.16 2.51007
8.28 0.383444
8.4 0.0336024
8.52 0.00138178
8.64 6.28082e-05
8.76 6.28082e-05
8.88 6.28082e-05
9 6.28082e-05
9.12 0
9.24 6.28082e-05
9.36 6.28082e-05
9.48 6.28082e-05
9.6 6.28082e-05
9.72 6.28082e-05
9.84 0
9.96 6.28082e-05
10.08 6.28082e-05
10.2 6.28082e-05
10.32 6.28082e-05
10.44 0
10.56 6.28082e-05
10.68 6.28082e-05
10.8 6.28082e-05
10.92 6.28082e-05
11.04 6.28082e-05
11.16 0
11.28 6.28082e-05
11.4 6.28082e-05
11.52 6.28082e-05
11.64 6.28082e-05
11.76 6.28082e-05
11.88 6.28082e-05
得到 $\D \g^{\text{lrc}}=0.0282 \; \text{kcal/mol-}\AA^2 = 19.592 \; \text{mJ/m}^2$
代码
| matlab | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | function 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
|
评论
- 2016-06-01 19:51:48
许小生x你好，我想问一下。如果我计算的是两相体系，求它们的界面张力，也是使用相同的方法吗？ - 2016-06-01 20:06:21
Jerkwin到下方的群里或论坛提问吧. - 2016-06-06 20:23:13
许小生x你这个测试用例的左侧是什么数据啊 - 2016-06-06 22:30:31
JerkwinTIP4P水的数密度分布, 第一列是z方向距离, nm, 第二列是数密度, 1/nm^3
