4.13. 隐式溶剂化计算¶
有时候，我们希望计算溶液中的体系性质，可能是溶液中的有机小分子，或者是插在电解质溶液中的电极。我们可以在构建晶胞的时候显式的把溶剂分子包括进去，但由于溶剂分子转动自由，影响需要平均化，因此需要包括非常大量的溶剂分子才能进行严格的计算，这样计算量太大了。
隐式溶剂化是一个半经验性的方法，把溶剂分子的影响抽象为一个介电常数的连续场，从而能以较小的计算量来考虑溶剂的影响，还可以通过简单的修正，考虑溶液中的电解质离子的影响。
目前，Hylanemos的隐式溶剂化还只能用于单点计算，不能用于结构优化。在后续的版本中，会增加结构优化的功能。
4.13.1. 二氧五环在水中的溶解热计算输入文件¶
计算二氧五环在水中的溶解热，可以分别计算真空中的二氧五环的能量，以及用隐式溶剂化计算水中的二氧五环的能量，相减得到。
真空中的计算是一个常规计算，这里不再讨论，只给出隐式溶剂化计算的输入JSON文件的参数设置。
JSON文件如下：
{
"job": {
"calculation_type": "scf",
"calc_force": true,
"calc_stress": true,
"occupation_method": "smearing",
"do_solvation": true
},
"job_io": {
"pp_dir": ".",
"pp_files": ["C_ONCV_PBE_sr.upf", "O_ONCV_PBE_sr.upf", "H_ONCV_PBE_sr.upf"],
"prefix": "solvation"
},
"cell": {
"lattice": [
10.0, 0.0, 0.0,
0.0, 10.0, 0.0,
0.0, 0.0, 10.0
],
"cell_units": "A"
},
"ions": {
"element_names": ["C", "O", "H"],
"element_nums": [3, 2, 6],
"positions": [
[0.3680194334, 0.6084306164, 0.5247952504],
[0.4545613772, 0.4919685329, 0.4799565398],
[0.5859967780, 0.6764796069, 0.5116746537],
[0.5849752609, 0.5343982457, 0.5214887499],
[0.4516843794, 0.7208641418, 0.4952392565],
[0.3459209776, 0.6029691618, 0.6327938762],
[0.2747083500, 0.6200258042, 0.4686169682],
[0.4325529933, 0.3970877785, 0.5299519883],
[0.4500386774, 0.4781772960, 0.3704700476],
[0.6435523285, 0.7106003098, 0.4238393414],
[0.6301294442, 0.7159785060, 0.6047233280]
]
},
"pw": {
"ecutwfc": 40.0
},
"kpts": {
"k_type": "automatic",
"k_mesh": [1, 1, 1]
},
"symmetry": {
"use_sym": false
},
"electron_step": {
"elec_max_steps": 100,
"elec_e_conv": 5e-7
},
"smearing": {
"smearing_alg": "gauss",
"smearing_width": 0.005
},
"solvation": {
"solvation_epsilon": 78,
"solvation_gamma": 50.0,
"solvation_P": -0.35
}
}
JSON输入参数介绍：
JSON输入文件可以分为几个模块，这里分别是job/job_io/pw/kpts/smearing/electron_step/cell/ions/symmetry/solvation， 大部分参数设置都和普通的结构优化计算相同，这里只介绍一些隐式溶剂化计算需要额外设置的参数。 下面将依次进行介绍。
job模块： 用来设置和本次计算类型相关的参数。
do_solvation：设置本次计算是否需要考虑隐式溶剂模型，这里设为true。
solvation模块： 用来设置和隐式溶剂化计算方法相关的参数。
solvation_epsilon：设置溶剂的介电常数，直接使用实验值，水可取78。solvation_gamma：设置溶剂的“表面张力”，但其实是若干形式相同的项的总和贡献，而不是单纯的表面张力，因此不能使用实验值，而是要对不同的体系拟合得到。水可取50.0。solvation_P：设置溶剂的“压强”，同样是若干形式相同的项的总和贡献，要对不同的体系拟合得到。水可取-0.35。（后面两个参数都是取自文献J. Chem. Phys. 136, 064102(2012)）
4.13.2. 执行计算¶
准备好JSON文件和赝势文件之后，按照 Hylanemos运行 中的方法执行计算。
4.13.3. 计算结果分析¶
在 relax_OUT 文件中的scf result模块可以看到计算结束后的体系能量。
真空中的计算结果为：
================================================================================
| scf result |
--------------------------------------------------------------------------------
Electron Step Finish
Total Energy: -54.62427262 Ha / -1486.40217908 eV
Fermi level: -1.50884866 eV
Energy by parts:
one_electron contribution = -166.04398597 Ha / -4518.28703107 eV
hartree contribution = 85.55594586 Ha / 2328.09588597 eV
xc contribution = -16.88958943 Ha / -459.58914094 eV
ewald contribution = 42.75335691 Ha / 1163.37810696 eV
smearing contribution = -0.00000000 Ha / -0.00000000 eV
隐式溶剂化计算结果为：
================================================================================
| scf result |
--------------------------------------------------------------------------------
Electron Step Finish
Total Energy: -54.62742482 Ha / -1486.48795486 eV
Fermi level: -1.38468181 eV
Energy by parts:
one_electron contribution = -166.14338513 Ha / -4520.99181992 eV
hartree contribution = 85.67911742 Ha / 2331.44755465 eV
xc contribution = -16.89815300 Ha / -459.82216756 eV
ewald contribution = 42.75335691 Ha / 1163.37810696 eV
smearing contribution = -0.00000000 Ha / -0.00000000 eV
solvation contribution = -0.01836103 Ha / -0.49962899 eV
因此溶解热为每个分子3.152e-3 Ha，换算得到1.97 kcal/mol。文献（J. Chem. Theory Comput. 2010, 6, 1509-1519）中给出实验值是4.10 kcal/mol，相差较大。
目前还没有实现带隐式溶剂化的结构优化，因此计算溶剂化效应的结构不是能量最低的结构，得到的溶解热偏小。后续增加这个功能后，能得到更好的结果。