计算化学公社

标题: cp2k运行AIMD设置NVT系综体系温度仍然持续升高问题 [打印本页]

作者
Author: zaq1 时间: 2024-12-26 14:16
标题: cp2k运行AIMD设置NVT系综体系温度仍然持续升高问题
本帖最后由 zaq1 于 2024-12-26 14:42 编辑

大家好，我用cp2k进行AIMD模拟。设置了NVT系综，但是看运行过程中体系的温度依旧持续升高，是因为我初始体系结构优化不合理吗，我试着解决这个问题，想着用退火，初始设置0.5fs，逐渐升高温度，达到300K就停止加温度，并持续到模拟结束，这个方法合理吗，但是cp2k好像不能设置控制体系温度的升高，需要什么命令大家知道嘛，关于设置的NVT系综如下，和输出的文件附上照片
&MOTION (, 下载次数 Times of downloads: 13) [url=]C3SOH.inp[/url]
  &MD
ENSEMBLE NVT
STEPS 100000 #Number of steps to run
TIMESTEP 1.0 #Step size in fs. Decrease it properly for high temperature simulation
TEMPERATURE 298.15 #Initial and maintained temperature (K)
# COMVEL_TOL 0 #Uncomment this can remove translation motion of center-of-mass every step
&THERMOSTAT
   TYPE NOSE
&END THERMOSTAT
&PRINT
   &PROGRAM_RUN_INFO
      &EACH
      MD    1 #Output frequency of MD information, 0 means never
      &END EACH
   &END PROGRAM_RUN_INFO
&END PRINT
  &END MD

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 14:35
分子动力学里热浴控温的原理就决定了体系的温度会围绕目标平衡值上下波动而不是一成不变，毕竟固定温度会给粒子动能分布带来过强的约束而造成不真实的运动。我在社长这帖的评论也说了，CP2K目前没有稳定速率线性升温的功能，温度设置的值就是目标的平衡温度。

温度上升幅度过大多少说明初始体系存在较大斥力，但是没有完整输入文件没法回答“初始体系结构优化不合理吗”的问题。另外Nose-Hoover热浴对尚未达到平衡的体系更是不推荐使用。

作者
Author: zaq1 时间: 2024-12-26 14:45

Uus/pMeC6H4-/キ发表于 2024-12-26 14:35
分子动力学里热浴控温的原理就决定了体系的温度会围绕目标平衡值上下波动而不是一成不变，毕竟固定温度会给 ...

好的老师，我上传了输入文件，但是这样的话温度就是一直在变化的了吗，cp2k不能保证体系温度一直不变是吗

作者
Author: zaq1 时间: 2024-12-26 14:49
本帖最后由 zaq1 于 2024-12-26 14:55 编辑

zaq1 发表于 2024-12-26 14:45
好的老师，我上传了输入文件，但是这样的话温度就是一直在变化的了吗，cp2k不能保证体系温度一直不变是吗

温度大概会在什么范围内波动呢，我这个温度上升的剧烈应该就是不合理的是吗，Nose-Hoover热浴不合适应该用什么，第二个问题，cp2k能否使用退火方法，如果能要用到什么命令呢

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 15:34
本帖最后由 Uus/pMeC6H4-/キ于 2024-12-26 15:42 编辑

zaq1 发表于 2024-12-26 14:45
好的老师，我上传了输入文件，但是这样的话温度就是一直在变化的了吗，cp2k不能保证体系温度一直不变是吗

不用叫老师哦……“温度一直在变化”有歧义，每步的瞬时温度涨落和一段时间的平均温度升降是两个不同的意思，参考下述社长博文。
谈谈怎么判断分子动力学模拟是否达到了平衡
http://sobereva.com/627（http://bbs.keinsci.com/thread-27122-1-1.html）

感觉基底的结构有问题，部分硅原子的[SiO4]四面体变形明显，靠近界面的钙原子配位没饱和。这是从某个体相直接切的slab么？

编辑：上面的帖子没刷出来，补充下回复。热浴改用CSVR就好，CP2K目前只能设置一个目标的平衡温度，在热浴作用下跑足够长时间动力学可以让平均温度符合这个值。

作者
Author: zaq1 时间: 2024-12-26 15:55
本帖最后由 zaq1 于 2024-12-26 15:57 编辑

Uus/pMeC6H4-/キ发表于 2024-12-26 15:34
不用叫老师哦……“温度一直在变化”有歧义，每步的瞬时温度涨落和一段时间的平均温度升降是两个不同的意 ...

好的，这个模型初始晶胞是斜的，不是正交，应该是我用了不正确的方法把它变为正交的才会变形，现在已改正（不知道是不是这个原因），还有就是温度的话，我设置为300K的话，模拟是随时间变化温度一直上升还是在300K附近波动，最后满足平均温度为目标值好呢，不知道什么样的变化属于正确的

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 17:03

zaq1 发表于 2024-12-26 15:55
好的，这个模型初始晶胞是斜的，不是正交，应该是我用了不正确的方法把它变为正交的才会变形，现在已改正 ...

有些斜的晶胞确实因为角度的数学问题不好改成正交的，是需要点注意。又仔细看了下输入文件，这计算如果用XYZ三维周期性的话基底相当于有上下两个表面，把基底移到盒子中间并在两面都覆盖上水层来与真空区隔开会更合理些，如果只要一个界面的话应该用XY二维周期性。水层可以摆得离基底近一些，毕竟表面有些钙原子因为配位不饱和会容易吸附水分子配位。

没有自定义的话，CP2K会利用温度设置值产生符合Maxwell分布的随机初速度，后续分子动力学模拟时热浴亦往此温度设置值调节体系粒子速度，瞬时温度的波动幅度由热浴的时间常数TIMECON控制。以我的经验，温度变化趋势和之前有没有对体系做几何优化（能量极小化）有关：
若做好再跑，瞬时温度会短暂下降（粒子偏离优化位置不多故受力较小、动能尚未提升到搭配热浴温度的程度）然后再逐渐上升向目标温度收敛；
而若不做就直接跑，瞬时温度确实会快速升高（粒子不合理接触产生作用力较大故动能猛增、热浴调节也没跟上）然后再逐渐下降向目标温度收敛。
考虑到后者容易让体系异常崩溃，还是推荐先做个几何优化减小受力再跑分子动力学模拟。

作者
Author: zaq1 时间: 2024-12-26 17:07
本帖最后由 zaq1 于 2024-12-26 17:09 编辑

Uus/pMeC6H4-/キ发表于 2024-12-26 17:03
有些斜的晶胞确实因为角度的数学问题不好改成正交的，是需要点注意。又仔细看了下输入文件，这计算如果用 ...

做了几何优化的，但可能由于前期那个结构正交化没做好导致结构优化能量不是最低，对了，请教一下您，怎么观察结构优化后得到的结构合理呢，收敛前提下，看键长键角吗，或者哪里有学习的资源吗，谢谢您，至于您说的做成上下面结构进行计算，我先把正交后的结构优化后先模拟一下观看一下温度的变化再考虑是否用这种结构，谢谢你啦

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 17:40

zaq1 发表于 2024-12-26 17:07
做了几何优化的，但可能由于前期那个结构正交化没做好导致结构优化能量不是最低，对了，请教一下您，怎么 ...

其实也没什么特别的，一方面注意SCF收敛精度，另一方面也结合化学常识判断分子取向、配位多面体、键长键角等几何特征。加水层之后跑分子动力学之前最好也做个几何优化，现在这个输入文件的水分子看起来比较稀疏随机，通过优化产生一些氢键可以让水层结构更像真实液相，对初期控温也有帮助。

我所知道的CP2K最佳的学习资源就是社长办的第一性原理计算培训班了。

作者
Author: zaq1 时间: 2024-12-26 19:05
本帖最后由 zaq1 于 2024-12-26 19:15 编辑

Uus/pMeC6H4-/キ发表于 2024-12-26 17:40
其实也没什么特别的，一方面注意SCF收敛精度，另一方面也结合化学常识判断分子取向、配位多面体、键长键 ...

好的，主要是不知道换了模型之后温度还会不会像这次跑的那么离谱，我先进行两次优化看一下，两个都是进行结构优化是嘛

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 19:56

zaq1 发表于 2024-12-26 19:05
好的，主要是不知道换了模型之后温度还会不会像这次跑的那么离谱，我先进行两次优化看一下，两个都是进行 ...

“两个”指什么？基底是实验测定的晶体结构截取的话，带真空层单独做一次结构优化能体现表面弛豫或表面重构的过程，在这个基础上再加水层做分子动力学模拟可能更贴近现实的干燥表面润湿，当然这也看具体研究课题。我上面说的分子动力学之前对水层和基底的整体优化，主要还是为了减小原子的受力免得模拟初期体系出现异常。

作者
Author: zaq1 时间: 2024-12-26 21:01
本帖最后由 zaq1 于 2024-12-26 21:26 编辑

Uus/pMeC6H4-/キ发表于 2024-12-26 19:56
“两个”指什么？基底是实验测定的晶体结构截取的话，带真空层单独做一次结构优化能体现表面弛豫或表面重 ...

我做的是AIMD模拟，就没做弛豫这部分，首先是对初始晶胞进行优化，结束后进行扩胞切面然后加模拟溶液，这个您说也可以优化，我想着也优化一次，然后再进行后续模拟，这里的初始晶胞指的是从晶体库下载的，这是重新结构优化得到的输出文件，您看是否有问题呢 (, 下载次数 Times of downloads: 6) [url=]C3S.out[/url]

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-26 22:22
本帖最后由 Uus/pMeC6H4-/キ于 2024-12-27 21:48 编辑

zaq1 发表于 2024-12-26 21:01
我做的是AIMD模拟，就没做弛豫这部分，首先是对初始晶胞进行优化，结束后进行扩胞切面然后加模拟溶液，这 ...

单看这个变胞优化输出文件，只能看出每步SCF收敛和总体位移、梯度（编辑：之前打错成受力了，更正一下）、压力收敛不错，看不出具体的结构变化如何。

晶体内部和表面的环境差异明显，直接把内部原子暴露出来不见得是稳定结构，表面弛豫与表面重构就是表面附近原子自发调整距离乃至改变成键方式适应新环境的过程。对当前体系而言这个因素很可能不能忽略，所以用变胞优化产生的结构切面扩胞添真空层之后，需要做一次几何优化让表面弛豫/重构稳定，之后添水层再做一次几何优化让水层结构更合理，这都会影响之后的分子动力学模拟。

作者
Author: zaq1 时间: 2024-12-27 09:19
本帖最后由 zaq1 于 2024-12-27 09:26 编辑

Uus/pMeC6H4-/キ发表于 2024-12-26 22:22
单看这个变胞优化输出文件，只能看出每步SCF收敛和总体受力、位移、压力收敛不错，看不出具体的结构变化 ...

所以无论做什么模拟，弛豫都是需要的是嘛，不是加了溶液进行模拟的，是变胞之后加真空层再进行优化，对了我还有一个问题，就是您说的之后再添水层做几何优化，是不论之后反应溶液加什么都要加个水层进行优化吗，还是直接加反应溶液后进行几何优化，因为我的反应溶液还加了OH根，不是纯水溶液，这是上面输出文件那个结构的轨迹文件，[url=] (, 下载次数 Times of downloads: 4) [/url]

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-27 10:43

zaq1 发表于 2024-12-27 09:19
所以无论做什么模拟，弛豫都是需要的是嘛，不是加了溶液进行模拟的，是变胞之后加真空层再进行优化，对了 ...

第一个问题，是。化学上的理由已经说了，毕竟表面弛豫这点很容易被问到，但是出自论坛上动画头像的陌生网友之口和出自现实中评委/审稿人之口肯定不是一个感觉吧（笑）

第二个问题，否。其他水溶液层和纯水层是一个操作流程，直接加到体系里和基底接触再优化、做动力学就行，体系该有啥就用啥。

注意CP2K产生的.xyz文件无晶胞参数信息，这样看变胞优化的轨迹只能检查盒子中间原子的成键情况，没法添加周期镜像验证边缘原子的跨盒子成键也合理。而扩胞也得确认用的晶胞参数是优化后的。

（另外帖子里有些多余的未解析url标签，用高级模式的纯文本编辑去掉为佳）

作者
Author: zaq1 时间: 2024-12-27 10:56

Uus/pMeC6H4-/キ发表于 2024-12-27 10:43
第一个问题，是。化学上的理由已经说了，毕竟表面弛豫这点很容易被问到，但是出自论坛上动画头像的陌生网 ...

哈哈，好的，还有最后一个小问题，就是我切面后加真空层再进行优化那输出得到的结构后面要怎么和我的溶液模型合并呢，已经有了真空层了，需要怎么操作呢，还有这几次优化的输入文件我都是用multiwfn软件生成的，任务类型都是选项4（优化结构和晶胞）对嘛，谢谢您啦

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-27 12:21

zaq1 发表于 2024-12-27 10:56
哈哈，好的，还有最后一个小问题，就是我切面后加真空层再进行优化那输出得到的结构后面要怎么和我的溶液 ...

真空层不影响建模啊，原来怎么加现在还怎么加呗

区分几何优化和变胞优化，几何优化（任务类型3）不改变晶胞参数只优化原子坐标，变胞优化（任务类型4）在特定外压下同时优化晶胞参数和原子坐标。有真空层的结构显然都该用几何优化，要不然外压直接把真空压没了。

作者
Author: zaq1 时间: 2024-12-27 13:37

Uus/pMeC6H4-/キ发表于 2024-12-27 12:21
真空层不影响建模啊，原来怎么加现在还怎么加呗

区分几何优化和变胞优化，几何优化（任务类型3）不改 ...

好的好的，谢谢你，我明白了

作者
Author: zaq1 时间: 2024-12-27 21:19
本帖最后由 zaq1 于 2024-12-27 21:23 编辑

Uus/pMeC6H4-/キ发表于 2024-12-27 12:21
真空层不影响建模啊，原来怎么加现在还怎么加呗

区分几何优化和变胞优化，几何优化（任务类型3）不改 ...

我用的是BFGS方法，我看到说用这个方法的话就必须有四个YES才算cp2k优化收敛，是这样的吗，我看我的都没有是还没收敛吗，我平常看是否收敛就是看SCF的步数，看Convergence的值，由这个方法判断是收敛了，可用YES这个方法就是不收敛的，这是为什么 (, 下载次数 Times of downloads: 5)

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-27 21:43
本帖最后由 Uus/pMeC6H4-/キ于 2024-12-27 21:49 编辑

zaq1 发表于 2024-12-27 21:19
我用的是BFGS方法，我看到说用这个方法的话就必须有四个YES才算cp2k优化收敛，是这样的吗，我看我的都没 ...

SCF迭代求解使波函数收敛（参考社长博文http://sobereva.com/665）和优化体系使位移、梯度等判据收敛是两码事，这张截图里对应后者。现在才刚结束第一步优化，体系离能量极小点的位置还远着呢，优化判据出现四个NO挺正常的，等后期在MAX_ITER定义的优化步数内能变成四个YES那才叫优化收敛。（编辑：四个YES对应的是几何优化收敛的情况，而变胞优化收敛还多一个压力的判据故应该有五个YES，12楼的输出文件就能看到。）

作者
Author: zaq1 时间: 2024-12-30 09:45
本帖最后由 zaq1 于 2024-12-30 10:50 编辑

Uus/pMeC6H4-/キ发表于 2024-12-27 21:43
SCF迭代求解使波函数收敛（参考社长博文http://sobereva.com/665）和优化体系使位移、梯度等判据收敛是两 ...

好的，我明白了，但是为什么对12楼那个体系进行扩胞切面加真空层进行几何优化了两天还没有收敛，我看了一下确实没达到4个YES，体系的总原子324个，是因为什么设置的不合理吗，C3S.zip[url=]C3S.inp[/url]

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-30 11:00

zaq1 发表于 2024-12-30 09:45
好的，我明白了，但是为什么对12楼那个体系进行扩胞切面加真空层进行几何优化了两天还没有收敛，我看了一 ...

和13楼一样，单看优化输出文件只能看出每步SCF收敛有点慢，看不出来具体输入的设置和输出的结构的变化情况。带比较大真空层的体系需要几百步才能几何优化收敛也不奇怪，如果结构还好的话继续跑下去看看。

作者
Author: zaq1 时间: 2024-12-30 11:10
本帖最后由 zaq1 于 2024-12-30 11:12 编辑

Uus/pMeC6H4-/キ发表于 2024-12-30 11:00
和13楼一样，单看优化输出文件只能看出每步SCF收敛有点慢，看不出来具体输入的设置和输出的结构的变化情 ...

那怎么可以看出输出结构的变化情况呢，能量现在是这个状态，跑的这种程度可以进行后面的计算了吗
(, 下载次数 Times of downloads: 6)

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-30 11:35

zaq1 发表于 2024-12-30 11:10
那怎么可以看出输出结构的变化情况呢，能量现在是这个状态，跑的这种程度可以进行后面的计算了吗

用可视化软件加载几何优化输出的-pos-1.xyz来观察呗，和后面做完分子动力学观察轨迹是一个道理。推荐使用1.9.3版的VMD，简单注册一下就能免费下载使用。

其实CP2K的优化并不直接以能量作为收敛判据，但是能量曲线后期这么平的话更应该看看结构究竟在如何变化了。

作者
Author: zaq1 时间: 2024-12-30 14:13

Uus/pMeC6H4-/キ发表于 2024-12-30 11:35
用可视化软件加载几何优化输出的-pos-1.xyz来观察呗，和后面做完分子动力学观察轨迹是一个道理。推荐使用 ...

观察轨迹文件就是看到原子的位置在轻微改变，原子坐标发生轻微变化，结构也没发生什么奇怪的变形，结构原子数324，带了15埃真空层，这个速度正常吗

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-30 15:49

zaq1 发表于 2024-12-30 14:13
观察轨迹文件就是看到原子的位置在轻微改变，原子坐标发生轻微变化，结构也没发生什么奇怪的变形，结构原 ...

刚看了下输入文件，才明白为什么计算速度慢了。现在的基底在z轴方向有点太厚了，即使是考虑表面弛豫/表面重构和后续分子动力学中水分子与基底的相互作用，本来也可以结合平移对称性选取等价的晶面构建薄一点的基底的。另外在Z方向加了真空层却还按XYZ的三维周期性来算的话，真空层的计算也会多花费一点时间，应该用XY二维周期性的，这样还可以把部分底层原子冻结进一步节约时间。

作者
Author: zaq1 时间: 2024-12-30 17:14
本帖最后由 zaq1 于 2024-12-30 17:15 编辑

Uus/pMeC6H4-/キ发表于 2024-12-30 15:49
刚看了下输入文件，才明白为什么计算速度慢了。现在的基底在z轴方向有点太厚了，即使是考虑表面弛豫/表面 ...

那我换一下周期性试一下，但是这里的基底厚度是我对初始斜晶胞进行计算正交后的厚度，我用的是MS软件建模，在MS里面Symmetry-redifine Lattice里面调整的，正交后自动就是这个厚度了，也不能改的吧应该，

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-30 18:26

zaq1 发表于 2024-12-30 17:14
那我换一下周期性试一下，但是这里的基底厚度是我对初始斜晶胞进行计算正交后的厚度，我用的是MS软件建模 ...

软件做正交化是要兼顾三维周期性的，但切表面加真空层以后只剩下二维周期性。当前基底明显可以按平移对称性划分出等价的三层（参考面心立方堆积的密排布层，我怀疑这个晶体本来就是三方晶系的），而且既然每层都不薄，已经超过可能的表面重构/表面弛豫的影响范围，取两层足够了。还可以设置下层原子冻结作为体相，只放开上层原子做几何优化和与水分子接触下的分子动力学，由于z方向一直有真空层故全用XY二维周期性。
(, 下载次数 Times of downloads: 8)
没有特别要求的话切表面可以考虑下原子的排布均衡，比如按划线的位置分隔层，让钙原子出现在环境等价但是周围有氧原子的位置，避免裸露钙原子的特殊电子结构使SCF难收敛的问题，这也有助于提高计算速度。

作者
Author: zaq1 时间: 2024-12-31 09:51
本帖最后由 zaq1 于 2024-12-31 11:32 编辑

Uus/pMeC6H4-/キ发表于 2024-12-30 18:26
软件做正交化是要兼顾三维周期性的，但切表面加真空层以后只剩下二维周期性。当前基底明显可以按平移对称 ...

1.切面之后我看是二维周期性了，但再加上真空层又是三维周期性了，后面也改不了
2.划分三层取两层的意思是就取两层的厚度加原子为模型再加真空层进行后续模拟吗（这个厚度的取法是否就是切面时将厚度调为两层的厚度再加真空层）
3.或者另一种方法，把最下面一层固定，还是两层固定只开放最上面一层模拟，后续进行AIMD最下面一层也是固定的，这个方法和2等效是吗

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-31 13:08

zaq1 发表于 2024-12-31 09:51
1.切面之后我看是二维周期性了，但再加上真空层又是三维周期性了，后面也改不了
2.划分三层取两层的意思 ...

我没用过你说的建模软件，不过把结构载入multiwfn之后可以先参考社长博文http://sobereva.com/610调整一下几何结构，再生成固定下层、按XY二维周期性计算的若干CP2K输入文件。从三层的结构删原子到保留二层也不难，这种操作用multiwfn、VMD乃至GaussView都能完成。

2和3的区别是指 “添加真空层+基底保留两层+固定下层放开上层” 和 “添加真空层+基底保留三层+固定下一/二层放开上层” 的区别么？我觉得不妨先用前者计算看看，毕竟计算资源有限，如果能用原子少一点的模型合理描述界面行为，没有必要用更大的。除非计算发现只放开一层不够，比如表面弛豫、表面重构影响较深，或者加溶液跑分子动力学时短时间内有明显的水分子渗透进入晶格、表面原子被剥落溶解进入溶液等现象，这样再考虑放开二层甚至更多。

作者
Author: zaq1 时间: 2024-12-31 13:27

Uus/pMeC6H4-/キ发表于 2024-12-31 13:08
我没用过你说的建模软件，不过把结构载入multiwfn之后可以先参考社长博文http://sobereva.com/610调整一 ...

好的我明白了，但是切面加了真空层的模型原本就是三维周期的，用multiwfn即使生成了二维周期计算的输入文件也会报错的，您的意思是用multiwfn把模型文件改为二维的吗，如果能改输入文件自然不是问题了

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-31 13:41

zaq1 发表于 2024-12-31 13:27
好的我明白了，但是切面加了真空层的模型原本就是三维周期的，用multiwfn即使生成了二维周期计算的输入文 ...

你是按http://sobereva.com/587进入创建CP2K输入文件的界面，然后用选项-7设定的XY周期性吗？具体是什么输入文件让CP2K报了什么错？

输入文件用文本编辑器是能修改，但鉴于其复杂性，不熟悉就最好别自己改动。我还是推荐报个CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）学一学。

作者
Author: zaq1 时间: 2024-12-31 18:39
本帖最后由 zaq1 于 2024-12-31 18:41 编辑

Uus/pMeC6H4-/キ发表于 2024-12-31 13:41
你是按http://sobereva.com/587进入创建CP2K输入文件的界面，然后用选项-7设定的XY周期性吗？具体是什么 ...

是的，按照-7设定的周期性，但是报错，报错显示如下，Index to radix array not found.

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-1-1 23:39

zaq1 发表于 2024-12-31 18:39
是的，按照-7设定的周期性，但是报错，报错显示如下，Index to radix array not found.

这个是因为基底太厚，用XY二维周期性时真空区也需要很大，z方向格点数量过多了。

那还是像原来那样用XYZ三维周期性好了，此时由于基底不是上下对称的，需要在-9 Other settings里用23 Set surface dipole correction开偶极校正。

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