cp2k跑硅酸盐熔体从头算 平衡结构不合理可能存在哪些原因

nusiew

大家好，我利用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基本都是游离的。我在想输入文件中是不是有一些参数设置的不合理才导致出现这种情况。我跑出来的平衡结构如下。

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


@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

sobereva

如果右图的结构的坐标是已知的，以其为初始结构，在低温下跑跑MD看看结构是否能维持。如果能，再上高温。
当前体系没必要加D3。把这种多余的设置去掉。

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

nusiew

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？

sobereva

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

删掉

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

nusiew

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

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

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

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

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

谢谢sob老师！

sobereva

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

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

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

2 NPT

759884279

你好，请问你的配位数相关问题解决了吗？我虽然是不同的熔盐体系，但也遇到了相关问题

nusiew

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

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

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

759884279

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

参考这个帖子 ...

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

759884279

你好，想再咨询一下你是如何统计总步数中阳离子与阴离子不同配位数占比的呢？

nusiew

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

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

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

759884279

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

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

nusiew

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

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

759884279

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

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