GROMACS表面张力单位,计算及其长程校正

类别:    标签: md gmx   阅读次数:   版权: (CC) BY-NC-SA

2014-09-24 13:00:48

单位

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

换算关系

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$, 则其表面张力的长程校正为

其中

$\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

评论

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 液滴接触角的计算
后一篇: 剪切模量各向异性

访问人次(2015年7月 9日起): | 最后更新: 2017-06-23 03:12:17 UTC | 版权所有 © 2008 - 2017 Jerkwin