计算化学公社

标题: 谈谈Dalton中的PCM溶剂模型的使用 [打印本页]

作者
Author: sobereva 时间: 2025-9-20 23:09
标题: 谈谈Dalton中的PCM溶剂模型的使用

谈谈Dalton中的PCM溶剂模型的使用

文/Sobereva@北京科音 2025-Sep-20

0 前言

Dalton程序有不可替代的价值，例如北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）里笔者会讲到使用Dalton算磷光速率、旋轨耦合矩阵元（包括激发态间的）、双光子吸收，这都是Dalton的特色功能。网上流传着一些错误显著错误的关于Dalton的溶剂模型的说法，甚至有的文中说什么“必须自定义空腔半径，否则计算无法继续下去”，明显误人子弟。本文专门说说用户应该了解的Dalton的PCM隐式溶剂模型的几个关键知识，其中不少内容在手册里都找不到。本文的情况对应Dalton 2018-2022版，对其它版本可能适用也可能不适用。

1 溶剂模型的使用

Dalton的隐式溶剂模型只支持PCM，不支持SMD等其它的。以下是个Dalton做B3LYP级别的单点+偶极矩计算的最简单的.dal文件例子，使用PCM模型描述水环境。

**DALTON INPUT
.RUN PROPERTIES
*PCM
.SOLVNT
WATER
*PCMCAV
**WAVE FUNCTIONS
.DFT
B3LYPg
**END OF DALTON INPUT

.SOLVNT下面是溶剂名，自带18种溶剂，包括（写化学式还是溶剂名都可以）：
H2O WATER
CH3OH METHANOL
C2H5OH ETHANOL
CHCL3 CLFORM
CH2CL2 METHYLCL
C2H4CL2 12DICLET
CCL4 TETRACLC
C6H6 BENZENE
C6H5CH3 TOLUENE
C6H5CL CLBENZ
CH3NO2 NITROMET
C7H16 N-EPTANE
C6H12 CYCLOHEX
C6H5NH2 ANILINE
CH3COCH3 ACETONE
THF TETHYDFU
DMSO DIMETSOX
CH3CN ACETONIT

输出文件中可以看到当前的溶剂信息
** LOOKING UP INTERNALLY STORED DATA FOR SOLVENT = WATER **
Optical and physical constants:
EPS= 78.390; EPSINF=  1.776; RSOLV=  1.385 A; VMOL=  18.070 ML/MOL;

*PCMCAV字段用于定义溶质孔洞，必须出现，如果下面没做额外设置，代表以默认的方式产生溶质的孔洞。Dalton用的是量子化学程序中普遍使用的靠原子范德华球叠加的方式构造孔洞。

*PCM里加入.NEQRSP可以使用非平衡溶剂，计算对应吸收的电子激发问题用得着。

如果要自定义新溶剂，比如静态介电常数为23.0、动态介电常数为1.84的溶剂，*PCM部分改为：
*PCM
.SOLVNT
WATER
.EPS
23
.EPSINF
1.84

2 Dalton构造溶质孔洞用的原子半径

Dalton带着PCM计算时，输出信息里可以看到诸如以下内容，给出的是各个原子的以Bohr为单位的X、Y、Z坐标和以埃为单位的构造溶质孔洞用的原子半径。

********SPHERES IN PCMSPHGEN************
INDEX       X       Y       Z       R
1 2.1499166192D+00 -1.9894550417D+00 1.6741801837D+00 1.2000000000D+00
2 4.0096560598D+00 2.5073910458D-01 0.0000000000D+00 1.2000000000D+00
3 2.1499166192D+00 -1.9894550417D+00 -1.6741801837D+00 1.2000000000D+00
...略

这些原子半径是怎么来的？实际上在源代码目录下的sirius/sircav.F中可以发现以下内容，可见用的是Bondi原子半径，但个C、N、O的半径是U. Pisa修改后的。对没定义半径的元素，半径直接当成0。
! A.Bondi, J.Phys.Chem. 68: 441-451(1964) gives alternate
! values, and a few transition metals.
      allocate(rvdw(99))
      rvdw = (/ 1.20d0, 1.22d0, 0.00d0, 0.00d0, 2.08d0, 1.85d0,
   &             1.54d0, 1.40d0, 1.35d0, 1.60d0, 2.31d0, 0.00d0,
   &             2.05d0, 2.00d0, 1.90d0, 1.85d0, 1.81d0, 1.91d0,
   &             2.31d0, 13*0.0d0, 2.00d0, 2.00d0, 1.95d0, 1.98d0,
   &             2.44d0, 13*0.0d0, 2.20d0, 2.20d0, 2.15d0, 0.00d0,
   &             2.62d0, 27*0.0d0, 2.40d0, 16*0.0d0 /)
! override the above table with U. Pisa's experience
! as to what works best for singly bonded C,N,O
      rvdw(6) = 1.70d0
      rvdw(7) = 1.60d0
      rvdw(8) = 1.50d0

所以，Dalton用的原子半径不是胡来的，是有依据的，不存在无条件的“必须自定义空腔半径”。仅当你的体系牵扯到部分元素没有内置半径时（而且那些元素的原子与溶剂有所接触时）才绝对需要补充半径。

3 Dalton中自定义半径实例

本节给一个例子，B3LYP/def2-SVP计算水中的乙醇，通过在Dalton中自定义原子半径使得其计算结果与Gaussian 16默认的IEFPCM的相吻合，从而加深读者对Dalton的PCM模型的理解和信心。

以下是Gaussian的输入文件，用scrf关键词默认的IEFPCM溶剂模型表现默认的水环境

#P B3LYP/def2SVP nosymm scrf
[空行]
Generated by Multiwfn
[空行]
  0  1
C    1.17229118 -0.41192328 0.00000000
H    1.13768688 -1.05277427 0.88593800
H    2.12181861 0.13268542 0.00000000
H    1.13768688 -1.05277427 -0.88593800
C    -0.00000000 0.55479430 0.00000000
H    0.05413938 1.20787442 -0.88657584
H    0.05413938 1.20787442 0.88657584
O    -1.19922622 -0.21255820 0.00000000
H    -1.94540849 0.40035379 0.00000000

输出文件在此：http://sobereva.com/attach/752/G16_scrf.out。从其中以下内容里可以看到溶剂信息和原子半径信息。其中Re0是原子半径（埃），Gaussian 16对IEFPCM模型默认用的是UFF力场定义的范德华半径。Alpha是给各个原子半径乘的系数。

Solvent             : Water, Eps=  78.355300 Eps(inf)= 1.777849
------------------------------------------------------------------------------
Spheres list:
ISph  on Nord    Re0 Alpha    Xe          Ye          Ze
1 C    1 1.9255  1.100    1.172291    -0.411923    0.000000
2 H    2 1.4430  1.100    1.137687    -1.052774    0.885938
3 H    3 1.4430  1.100    2.121819    0.132685    0.000000
4 H    4 1.4430  1.100    1.137687    -1.052774    -0.885938
5 C    5 1.9255  1.100    0.000000    0.554794    0.000000
6 H    6 1.4430  1.100    0.054139    1.207874    -0.886576
7 H    7 1.4430  1.100    0.054139    1.207874    0.886576
8 O    8 1.7500  1.100    -1.199226    -0.212558    0.000000
9 H    9 1.4430  1.100    -1.945408    0.400354    0.000000
------------------------------------------------------------------------------

当前Gaussian 16计算的偶极矩为1.8592 D。

下面是Dalton的.dal文件。.ICESPH设2代表要自定义原子半径，.NESFP设9代表要读入9个原子的半径（乙醇有9个原子，此例全都自定义）。.INA下面写上各个要修改半径的原子的序号。.RIN下面是各个原子的半径（埃），这里改成和Gaussian用的相同的UFF半径。注意顺序要按照.mol文件里的写，本例的.mol见http://sobereva.com/attach/752/ethanol.mol。.ALPHA定义各个原子的alpha值，都设成和Gaussian相同的1.1。Dalton里的B3LYP关键词对应的并不是Gaussian里的B3LYP，因此这里写B3LYPg使用和Gaussian相同的定义。

**DALTON INPUT
.RUN PROPERTIES
*PCM
.SOLVNT
WATER
.ICESPH
2
.NESFP
9
*PCMCAV
.INA
1
2
3
4
5
6
7
8
9
.RIN
1.4430
1.4430
1.4430
1.4430
1.4430
1.4430
1.9255
1.9255
1.7500
.ALPHA
1.1
1.1
1.1
1.1
1.1
1.1
1.1
1.1
1.1
**WAVE FUNCTIONS
.DFT
B3LYPg
**END OF DALTON INPUT

Dalton的输出文件是http://sobereva.com/attach/752/PCMmodcav_ethanol.out，可见其中的偶极矩的计算结果为1.8599 D，和Gaussian 16给出的1.8592完美吻合！

再顺带一提怎么让Gaussian得到Dalton默认孔洞设置下的结果。Dalton对H、C、O设的半径分别是1.2、1.7、1.5埃，alpha默认用的是1.2，因此Gaussian输入文件要写成下面这样与之对应：

#P B3LYP/def2SVP nosymm scrf=read
[空行]
Generated by Multiwfn
[空行]
  0  1
C    1.17229118 -0.41192328 0.00000000
H    1.13768688 -1.05277427 0.88593800
H    2.12181861 0.13268542 0.00000000
H    1.13768688 -1.05277427 -0.88593800
C    -0.00000000 0.55479430 0.00000000
H    0.05413938 1.20787442 -0.88657584
H    0.05413938 1.20787442 0.88657584
O    -1.19922622 -0.21255820 0.00000000
H    -1.94540849 0.40035379 0.00000000
[空行]
alpha=1.2
modifysph
[空行]
H 1.2
C 1.7
O 1.5
[空行]
[空行]

Gaussian 16的偶极矩计算结果为1.9136 D，Dalton在默认孔洞下计算的输出文件为http://sobereva.com/attach/752/PCM_ethanol.out，结果为1.9139 D，可见再次吻合得极好！

4 关于Dalton用PCM计算时的警告

直接按照前面说的方式在Dalton中用PCM模型做计算，最终会出现警告

576: Warning, element          34       263 of SI too big: set to zero
577: Warning, element          50       249 of SI too big: set to zero
578: Warning, element          52       249 of SI too big: set to zero
579: Warning, element          56       247 of SI too big: set to zero
...

很多人看到这警告就大惊失色了。实际上这些Warning可以直接无视，如果想不显示这些Warning，在*PCM里加上以下内容即可（在Dalton自带的涉及PCM的测试文件里总是带着这个）
.NPCMMT
0

这里说一下原因。.NPCMMT选项在手册里找不到，但在源代码sirius/sirpcm.F里能看到其说明
NPCMMT = 0 No correction of the DI, SI and C matrices
NPCMMT = 1 Correction of DI and SI  (default)
NPCMMT = 2 Correction of DI, SI and C
进而看sirief.F，会发现NPCMMT为1、2的时候都会自动做矩阵的检查，并可能导致那些warning的出现。是否会出现warning和是否自定义原子半径并没必然关系。NPCMMT为1、2时做的校正其实没什么意义，和为0时的结果差异微乎其微。对于前例，为1时偶极矩为1.913867、单点能为-154.9289657907 Ha，为0时分别为1.913182和-154.9289638876 Ha。

计算时还有可能出现** WARNING ** A VERY POOR TESSELATION HAS BEEN CHOSEN，这通过设NPCMMT为0也不会避免。根据我的测试，这个warning一般也是无害的，结果还是正常的。实在不放心的话可以按照上文的方式，把Gaussian中的原子半径设成与Dalton相同后对Dalton的结果做验证，也可以尝试把Dalton的原子半径设成和Gaussian默认的一致看看是否warning能消除，且与Gaussian结果相符。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)