What are the possible reasons for the producted unreasonable structure of silicate melt using cp2k ab initio MD?

[Sun Yiwei]

Hello everyone, I use cp2k to run an ab initio calculation of the silicate melt structure. I tried many possible situations, but the structure that came out was unreasonable. I don't know where the problem occurred. Please help me to see the possible reasons. The inp file is as follows.

The system contains 32 Ca, 50 Si, 138 O, 4 B, a total of 224 atoms, and the atoms are randomly distributed in the initial structure. Theoretically speaking, there should be 4 O atoms coordinated around Si atoms in the equilibrium structure, but in the structure I ran out, only part of the Si atoms are surrounded by 4 O atoms, and the other part of the Si atoms has only three O atoms, two O atoms, or one O atom coordinated, and even a few Si atoms are free. In addition, in the structure I ran out, the Si inside the system was basically 4-coordinated, and the coordination numbers of Si and O atoms at the system boundary did not reach saturation, as shown in the figure I simulated below. In theory, O atoms should be combined with Si atoms, but the O atoms at the edges of the system are basically free. I was wondering if there are some unreasonable parameter settings in inp file that caused this situation. The equilibrium structure I ran out is as follows.

I have been stuck here for a long time, so please help. Thank you all for your help!

@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 800        

            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

            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 10000

        TEMPERATURE 1823

        TIMESTEP 2.0

    &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

[Thomas Kühne]

Dear Sun Yiwei, 

in my opinion your temperature is too low, possibly even below the experimental melting temperature. Here, you have to be aware 

that the temperature scale is rather fine and DFT easily mispredicts the transition temperature by ~200 K. Moreover, close to the 

transition line the time-scales for melting/crystallization become unreasonably high for direct MD simulations. Hence, I conjecture 

that what you obtained may not a genuine melt and the eventual structure strongly dependent on your chosen random initial 

structure. Also, be aware that for pure silica, the particle density of the liquid is approx. 20% below that of its crystalline phase. 

Instead, I would suggest to equilibrate the melt at a much higher temperature and possibly then quench it to the desired target 

temperature well above the experimental transition temperature …

Best, 

Thomas Kühne

Am 07.09.2021 um 09:56 schrieb Sun Yiwei <sun789...@gmail.com>:

Hello everyone, I use cp2k to run an ab initio calculation of the silicate melt structure. I tried many possible situations, but the structure that came out was unreasonable. I don't know where the problem occurred. Please help me to see the possible reasons. The inp file is as follows.

The system contains 32 Ca, 50 Si, 138 O, 4 B, a total of 224 atoms, and the atoms are randomly distributed in the initial structure. Theoretically speaking, there should be 4 O atoms coordinated around Si atoms in the equilibrium structure, but in the structure I ran out, only part of the Si atoms are surrounded by 4 O atoms, and the other part of the Si atoms has only three O atoms, two O atoms, or one O atom coordinated, and even a few Si atoms are free. In addition, in the structure I ran out, the Si inside the system was basically 4-coordinated, and the coordination numbers of Si and O atoms at the system boundary did not reach saturation, as shown in the figure I simulated below. In theory, O atoms should be combined with Si atoms, but the O atoms at the edges of the system are basically free. I was wondering if there are some unreasonable parameter settings in inp file that caused this situation. The equilibrium structure I ran out is as follows.

I have been stuck here for a long time, so please help. Thank you all for your help!

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/dad379b2-a39a-413c-96c0-507d82db1b71n%40googlegroups.com.
<1555.png>