计算化学公社

标题: cp2k跑硅酸盐熔体从头算平衡结构不合理可能存在哪些原因 [打印本页]

作者
Author: nusiew 时间: 2021-9-7 13:12
标题: cp2k跑硅酸盐熔体从头算平衡结构不合理可能存在哪些原因
本帖最后由 nusiew 于 2021-9-12 15:57 编辑

大家好，我利用cp2k跑硅酸盐熔体结构的从头算，做了很多尝试，但跑出来的结构一直不合理，不知道究竟哪里出现了问题，请大家帮我看一看可能出现的原因，输入的控制文件如下。

体系中包含32个Ca, 50个Si, 138个O, 4个B，共224个原子，初始结构中原子随机均匀分布。理论上得到的平衡结构中Si的周围应该有4个O配位，但我跑出来的结构中只有一部分Si的周围是4个O配位，另外一部分Si的配位包含3配位，2配位，1配位，甚至少数Si是游离的。除此之外，我跑出来的结构中，内部的Si基本都是4配位，外部的Si和O配位数都是不饱和的，理论上O应该与Si结合，但是外部的O基本都是游离的。我在想输入文件中是不是有一些参数设置的不合理才导致出现这种情况。我跑出来的平衡结构如下。

因为在这里已经卡了很久了，一直没有突破，请大家帮忙看看。抱拳了，各位老铁！

(, 下载次数 Times of downloads: 82)
@SET SYSNAME CaSiOB
@SET CELL 14.7599

&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT_NAME ${SYSNAME}
RUN_TYPE MD
&END GLOBAL

&FORCE_EVAL
METHOD  QS
STRESS_TENSOR  ANALYTICAL
&DFT
      BASIS_SET_FILE_NAME BASIS_MOLOPT
      POTENTIAL_FILE_NAME GTH_POTENTIALS
      &SCF
         MAX_SCF 200
         EPS_SCF 1.0e-5
         SCF_GUESS  ATOMIC
         &PRINT
            &RESTART OFF
            &END RESTART
         &END PRINT
         &OT ON
            PRECONDITIONER  FULL_SINGLE_INVERSE
            MINIMIZER    DIIS
         &END OT
         &OUTER_SCF
            EPS_SCF 1.0e-5
            MAX_SCF  100
         &END OUTER_SCF
      &END SCF
      &QS
         EPS_DEFAULT    1e-10
      &END QS
      &MGRID
         NGRIDS 4
         CUTOFF 500
         REL_CUTOFF 60
      &END MGRID
      &XC
         &XC_FUNCTIONAL PBE
         &END XC_FUNCTIONAL
         DENSITY_CUTOFF    1e-10
         GRADIENT_CUTOFF    1e-10
         TAU_CUTOFF       1e-10
         &VDW_POTENTIAL
            POTENTIAL_TYPE  PAIR_POTENTIAL
            &PAIR_POTENTIAL
                  R_CUTOFF 18
                  TYPE       DFTD3
                  PARAMETER_FILE_NAME dftd3.dat
                  REFERENCE_FUNCTIONAL PBE
            &END PAIR_POTENTIAL
         &END VDW_POTENTIAL
      &END XC
&END DFT
&SUBSYS
      &TOPOLOGY
         COORD_FILE ${SYSNAME}.xyz
         COORD_FILE_FORMAT XYZ
         CONN_FILE_FORMAT OFF
      &END TOPOLOGY
      &CELL
         &CELL_REF
            ABC ${CELL} ${CELL} ${CELL}
            ALPHA_BETA_GAMMA 90.0 90.0 90.0
            MULTIPLE_UNIT_CELL  1 1 1
            PERIODIC XYZ
         &END CELL_REF
         ABC ${CELL} ${CELL} ${CELL}
         ALPHA_BETA_GAMMA 90.0 90.0 90.0
         MULTIPLE_UNIT_CELL  1 1 1
         PERIODIC XYZ
      &END CELL
      &KIND B
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE-q3
      &END KIND
      &KIND O
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE-q6
      &END KIND
      &KIND Si
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE-q4
      &END KIND
      &KIND Ca
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE-q10
      &END KIND
&END SUBSYS
&END FORCE_EVAL

&MOTION
&MD
      ENSEMBLE NVT
      &THERMOSTAT
         REGION MASSIVE
         TYPE NOSE
         &NOSE
            TIMECON 100
         &END NOSE
      &END THERMOSTAT
      STEPS 20000
      TEMPERATURE 2500
      TIMESTEP 0.5
&END MD
&PRINT
      &RESTART
         FILENAME =${SYSNAME}.rst
         &EACH
            MD 10
         &END EACH
      &END RESTART
      &STRESS
         FILENAME =${SYSNAME}.stress
         &EACH
            MD 10
         &END EACH
      &END STRESS
      &TRAJECTORY
         FILENAME =${SYSNAME}.pos.xyz
         FORMAT XYZ
         UNIT angstrom
         &EACH
            MD 10
         &END EACH
      &END TRAJECTORY
      &VELOCITIES
         FILENAME =${SYSNAME}.vel.xyz
         FORMAT XYZ
         UNIT angstrom/fs
         &EACH
            MD 10
         &END EACH
      &END VELOCITIES
      &FORCES
         FILENAME =${SYSNAME}.frc.xyz
         FORMAT XYZ
         UNIT amu*angstrom/fs^2
         &EACH
            MD 10
         &END EACH
      &END FORCES
      &CELL
         FILENAME =${SYSNAME}.cell
         &EACH
            MD 10
         &END EACH
      &END CELL
&END PRINT
&END MOTION

作者
Author: sobereva 时间: 2021-9-8 04:03
如果右图的结构的坐标是已知的，以其为初始结构，在低温下跑跑MD看看结构是否能维持。如果能，再上高温。
当前体系没必要加D3。把这种多余的设置去掉。

体系达到平衡之前不建议用Nose-Hoover热浴。建议用CSVR

作者
Author: nusiew 时间: 2021-9-8 11:25
本帖最后由 nusiew 于 2021-9-8 17:14 编辑

sobereva 发表于 2021-9-8 04:03
如果右图的结构的坐标是已知的，以其为初始结构，在低温下跑跑MD看看结构是否能维持。如果能，再上高温。
...

当前体系没必要加D3是不是意味着下面的范德华校正部分可以删掉？
         &VDW_POTENTIAL
            POTENTIAL_TYPE  PAIR_POTENTIAL
            &PAIR_POTENTIAL
                  R_CUTOFF 18
                  TYPE       DFTD3
                  PARAMETER_FILE_NAME dftd3.dat
                  REFERENCE_FUNCTIONAL PBE
            &END PAIR_POTENTIAL
         &END VDW_POTENTIAL

我用经典分子动力学跑出来的平衡结构作为初始结构来跑从头算，在2500K条件下很快就能收敛且结构基本没有太大变化。由于缺少某些元素力场，目前我是想直接用cp2k基于原子随机分布的初始结构中跑出合理的平衡结构，但一直跑不出来。由于目标温度在1823K，初始结构结构是混乱的，是不是应该先在高于1823K的温度下跑然后再逐步降温到1823K？

作者
Author: sobereva 时间: 2021-9-8 23:24

nusiew 发表于 2021-9-8 11:25
当前体系没必要加D3是不是意味着下面的范德华校正部分可以删掉？
&VDW_POTENTIAL[/backcolo ...

删掉

可以试试更高温度然后逐渐降温

作者
Author: nusiew 时间: 2021-9-9 10:42

sobereva 发表于 2021-9-8 23:24
删掉

可以试试更高温度然后逐渐降温

谢谢sob老师！还有一个些小疑问就是：

1. 我这个体系为啥没必要加范德华校正呢？是因为范德华校正对原子间作用力贡献太小了吗？

2. 因为我最开始是参考1823K的熔体密度建立的随机均匀分布初始结构模型，考虑到不同温度下熔体密度的变化，如果我从更高温如3500K开始逐渐降温的话，那么更高温以及降温过程选用什么系综比较合适？是用NVT还是NPT呢？

谢谢sob老师！

作者
Author: sobereva 时间: 2021-9-10 02:21

nusiew 发表于 2021-9-9 10:42
谢谢sob老师！还有一个些小疑问就是：

1. 我这个体系为啥没必要加范德华校正呢？是因为范德华校正对原 ...

1 对此体系考虑色散校正的必要性很低。而且当前遇到莫名其妙的问题，自然尽可能先把计算设置简化便于找原原因

2 NPT

作者
Author: 759884279 时间: 2024-4-23 09:15
你好，请问你的配位数相关问题解决了吗？我虽然是不同的熔盐体系，但也遇到了相关问题

作者
Author: nusiew 时间: 2024-4-23 21:52

759884279 发表于 2024-4-23 09:15
你好，请问你的配位数相关问题解决了吗？我虽然是不同的熔盐体系，但也遇到了相关问题

cp2k分子动力学输出的是原子的真实轨迹，需要进行周期性变换，用VMD很容易就转换过来了

参考这个帖子：http://bbs.keinsci.com/thread-24834-1-1.html

作者
Author: 759884279 时间: 2024-5-9 11:14

nusiew 发表于 2024-4-23 21:52
cp2k分子动力学输出的是原子的真实轨迹，需要进行周期性变换，用VMD很容易就转换过来了

参考这个帖子 ...

感谢你的回复，之前因为等级不够一直无法回复。我也用VMD进行了分析，但我的是配比总是大于已有文献报道，不知道该如何是好

作者
Author: 759884279 时间: 2024-5-13 08:36
你好，想再咨询一下你是如何统计总步数中阳离子与阴离子不同配位数占比的呢？

作者
Author: nusiew 时间: 2024-5-20 16:07

759884279 发表于 2024-5-13 08:36
你好，想再咨询一下你是如何统计总步数中阳离子与阴离子不同配位数占比的呢？

给不同的阴-阳离子对设置不同的最近邻截断半径，使用ISAACS、RINGS等软件，或者自己写脚本即可计算出来

https://sourceforge.net/projects/isaacs/
https://sourceforge.net/projects/rings-code/

作者
Author: 759884279 时间: 2024-5-20 22:07

nusiew 发表于 2024-5-20 16:07
给不同的阴-阳离子对设置不同的最近邻截断半径，使用ISAACS、RINGS等软件，或者自己写脚本即可计算出来
...

感谢你的回复，RINGS这款软件CP2K的输出轨迹文件.xyz也可以使用吗？之前用VASP的同学倒是提到过，没想到CP2K也可以，受教了

作者
Author: nusiew 时间: 2024-5-21 08:55

759884279 发表于 2024-5-20 22:07
感谢你的回复，RINGS这款软件CP2K的输出轨迹文件.xyz也可以使用吗？之前用VASP的同学倒是提到过，没想到C ...

xyz轨迹文件的格式通常是通用的，而且比较简单

作者
Author: 759884279 时间: 2024-5-21 09:34

nusiew 发表于 2024-5-21 08:55
xyz轨迹文件的格式通常是通用的，而且比较简单

感谢，最后想在问问论坛里或者网上有RINGS这款软件的相关教程吗？从零碎的消息中看到似乎与社长开发的Multiwfn类似

作者
Author: Jiangqin 时间: 2025-6-24 15:51
老师，您好！请问您右图的理论xyz文件是怎么生成的哦，是用分子动力学软件得到的吗，不用软件可以直接生成吗？

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