GROMACS表面张力单位,计算及其长程校正 |
您所在的位置：网站首页 › mn和kg换算 › GROMACS表面张力单位,计算及其长程校正 |
|
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 47function 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(ksiRc)); 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 Jerkwin TIP4P水的数密度分布, 第一列是z方向距离, nm, 第二列是数密度, 1/nm^3 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |