Modeling stishovite with ReaxFF

LAMMPS users, I am able to reproduce the density as reported in 2003 ReaxFF paper for quartz, coesite, and cristobalite. Yet, for stishovite, I am not able to do so. Initially, I use NPT to equilibrate at 1 atm and 300 K, like I did with other polymorphs.

But with stishovite, the temperature shoots up to ~2,000 K before going back to 300 K. However, the density drops to about ~3.2 g/cc, much lower than what Dr. van Duin reported of ~4.29 g/cc.

I know my structure is stable because when I use other potentials (Tersoff, BKS, COMB10), the structure is stable, and by that, I mean the temperature does not shoot up, and the density calculated is reasonable compared to exp.

I am using a timestep of 0.2 fs upon the suggesting of Dr. van Duin, and even when changing this to 0.1 fs, or even 0.05 fs, I get the same results.

I am using real units.

My npt ensemble looks like:

fix 1 all npt temp 300 300 20 x 1 1 200 y 1 1 200 z 1 1 200

Does anyone know what could be the culprit behind the hike in temp from 300 K to 2000 K, before it drops slowly back down to 300 K?

Note I am also modeling ~768 atoms, but no change when modeling ~3,000 atoms.

Thanks,

Ben

LAMMPS users, I am able to reproduce the density as reported in 2003 ReaxFF
paper for quartz, coesite, and cristobalite. Yet, for stishovite, I am not
able to do so. Initially, I use NPT to equilibrate at 1 atm and 300 K, like
I did with other polymorphs.

But with stishovite, the temperature shoots up to ~2,000 K before going back

which means, that your geometry is not commensurate with the force
field parameters. your system is at very high potential energy and
"wants" to go to the next minimum. i would very carefully, double,
triple and quadruple check, if the force field parameters and all
other settings (units!) are correct.

axel.

That is exactly what I have done. I always check these parameters and geometry.

I triple checked the geometry by running the same geometry with different FFs and had no issues.
I triple checked the FF by running the same ffield.reax file with other geometries (quartz, coesite, cristobalite), and had no issues.

So, by doing this, I (feel) like I have shown that the geometry is not the issue and the FF is not the issue.

I know your advice makes sense, and the other 99 times I have had issues, your advice regarding just checking geom and FF has worked. For some reason, this time, specifically for this polymorph, and this potential, I cannot figure out what is different.

Ben

Is there a track record of ReaxFF SiO being able to reproduce correctly dynamics of stishovite? In the original paper (van Duin, J phys chem A, 2003), it looks like all that has been done for stishovite (and some other polymorphs) are single point energy calculations at different densities, to check the equation of state. I don't want to criticize ReaxFF, which is a brilliant approach to reactive forcefields, but this is a typical flaw often present in papers introducing new parametrizations of forcefields. It is not enough to compute the energy of a system, keeping its structural aspect constant and only varying its density (isotropically in general). It's not enough either to only perform energy minimization. The only correct way is to perform dynamics at a specific pressure and temperature.

I've tried to reproduce your results using reax/c, and using the FF parameters file included in the supporting information of the J phys chem C paper (Kulkarni, 2012). It's important to note that this is a reparametrization of the original ReaxFF SiO forcefield, which is also a possible explanation of why stishovite isn't correctly described. Anyway, I'm getting something very similar to what you get:

Step time Press Volume sysdensi Temp PotEng KinEng Lx
        0 1e-06 2376606.1 11797.19 4.2286509 300 -90883.09 1340.4703 20.967055
        5 0.500001 2295122.7 11998.305 4.1577703 316.08153 -97816.608 1412.3263 21.085531
       10 1.000001 2094389.3 12496.894 3.9918877 359.0828 -113785.46 1604.4661 21.373646
       15 1.500001 3092071.1 12991.92 3.8397864 496.82117 -130738.42 2219.9134 21.652216
       20 2.000001 2602151.3 13649.352 3.6548399 850.6331 -168802.81 3800.828 22.011447
       25 2.500001 1108334.2 14224.98 3.5069432 653.24188 -182432.01 2918.8377 22.316623
       30 3.000001 799752.45 14596.52 3.4176775 718.05464 -187777.57 3208.4363 22.50925
       35 3.500001 697018.21 14886.064 3.3512014 788.16818 -191277.27 3521.7201 22.657111
       40 4.000001 618450.93 15144.747 3.2939605 873.19227 -194205.9 3901.6276 22.7876
       45 4.500001 553249.73 15388 3.2418897 949.65848 -196693.37 4243.2966 22.908956
       50 5.000001 506973.13 15620.414 3.193654 1017.0626 -198873.61 4544.4738 23.023716
       55 5.500001 472983.82 15845.309 3.148326 1081.8947 -200873.67 4834.159 23.133685
       60 6.000001 455643.65 16065.651 3.1051463 1147.8889 -202778.67 5129.0366 23.240423
       65 6.500001 465353.55 16285.645 3.0632007 1212.6409 -204676.22 5418.3635 23.346022
       70 7.000001 508999.15 16512.994 3.0210269 1282.5126 -206753.57 5730.567 23.454158
       75 7.500001 556936.19 16758.371 2.9767927 1366.3017 -209229.15 6104.9562 23.569761
       80 8.000001 517261.19 17024.543 2.9302517 1439.5119 -211884.08 6432.0766 23.693892
       85 8.500001 356472.93 17289.812 2.8852943 1490.9041 -214071.38 6661.709 23.816321
       90 9.000001 125619.9 17513.817 2.848391 1521.8579 -215257.22 6800.0178 23.918733
       95 9.500001 -66047.541 17661.942 2.8245024 1539.7145 -215669.06 6879.8051 23.985976
      100 10.000001 -180872.43 17722.762 2.8148095 1550.5969 -215881.94 6928.4304 24.013477
      105 10.500001 -204381.6 17705.206 2.8176005 1561.3951 -216261.67 6976.679 24.005545
      110 11.000001 -161765.41 17634.636 2.828876 1579.3993 -216821.31 7057.1262 23.973608
      115 11.500001 -91567.102 17543.17 2.843625 1608.6962 -217422.13 7188.0314 23.932089
      120 12.000001 -31383.975 17458.075 2.8574856 1648.8249 -217992 7367.3357 23.893331
      125 12.500001 -2126.1199 17393.485 2.8680967 1695.8716 -218549.39 7577.5514 23.863828
      130 13.000001 -6422.7306 17349.718 2.8753318 1744.2579 -219119.32 7793.7528 23.843795
      135 13.500001 -35987.454 17317.638 2.8806583 1788.0754 -219696.58 7989.54 23.82909
      140 14.000001 -73520.572 17284.223 2.8862274 1822.9443 -220270.25 8145.342 23.813754
      145 14.500001 -99156.844 17238.572 2.8938707 1848.1032 -220848.25 8257.758 23.79277
      150 15.000001 -102349.11 17176.711 2.9042927 1866.3617 -221442.3 8339.3413 23.764276
      155 15.500001 -86921.071 17102.464 2.9169011 1881.6086 -222045.3 8407.4682 23.729986
      160 16.000001 -65287.859 17023.992 2.9303466 1895.8234 -222638.04 8470.9831 23.693636
      165 16.500001 -50234.468 16948.772 2.9433517 1908.0169 -223207.48 8525.4665 23.658688
      170 17.000001 -46804.922 16880.123 2.9553219 1915.3091 -223749.33 8558.0499 23.626702
      175 17.500001 -53158.533 16817.009 2.9664132 1914.5959 -224260.08 8554.863 23.597219
      180 18.000001 -64105.968 16755.899 2.9772318 1904.0306 -224737.39 8507.6547 23.568602
      185 18.500001 -72932.372 16693.056 2.9884401 1883.4485 -225181.16 8415.6891 23.5391
      190 19.000001 -74170.968 16626.518 3.0003995 1854.0341 -225590.88 8284.2586 23.507783
      195 19.500001 -67485.428 16557.058 3.0129869 1817.399 -225961.81 8120.5647 23.475002
      200 20.000001 -56677.094 16487.442 3.0257087 1774.6076 -226286.76 7929.3625 23.442054
      205 20.500001 -46294.67 16420.785 3.037991 1725.9138 -226560.05 7711.7873 23.410421
      210 21.000001 -40862.41 16359.097 3.0494468 1671.2424 -226781.24 7467.5025 23.381069
      215 21.500001 -41248.698 16302.442 3.0600444 1610.6708 -226954.2 7196.8547 23.354046
      220 22.000001 -44683.249 16249.297 3.0700527 1544.8357 -227084.44 6902.6877 23.328641
      225 22.500001 -47078.827 16197.747 3.0798232 1475.2262 -227178.71 6591.6563 23.303945
      230 23.000001 -46221.037 16146.678 3.0895641 1403.9345 -227244.09 6273.1084 23.279428
      235 23.500001 -42174.478 16096.157 3.0992613 1332.632 -227285.13 5954.5119 23.255123
      240 24.000001 -36342.59 16047.075 3.1087408 1262.3834 -227304.41 5640.6246 23.231462
      245 24.500001 -31807.269 16000.509 3.1177881 1193.7434 -227304.37 5333.9254 23.208969
      250 25.000001 -30632.719 15956.907 3.1263074 1126.8988 -227288.26 5035.2479 23.187868
      255 25.500001 -32324.044 15915.743 3.1343933 1062.176 -227260.83 4746.0512 23.167911
      260 26.000001 -34390.484 15875.907 3.1422582 1000.2214 -227227.37 4469.2234 23.148566
      265 26.500001 -33895.458 15836.513 3.1500746 941.60868 -227191.31 4207.3282 23.129403
      270 27.000001 -29536.789 15797.553 3.1578434 886.67943 -227153.51 3961.8914 23.11042
      275 27.500001 -22145.132 15760.027 3.1653625 835.59822 -227113.55 3733.6486 23.092107
      280 28.000001 -14745.987 15725.443 3.1723239 788.28942 -227071.39 3522.2618 23.075203
      285 28.500001 -10660.908 15694.902 3.1784968 744.57596 -227028.73 3326.9398 23.060256
      290 29.000001 -10080.629 15668.444 3.1838642 704.48838 -226988.96 3147.8191 23.04729
      295 29.500001 -10793.703 15645.25 3.1885843 668.31088 -226955.7 2986.1696 23.035912
      300 30.000001 -9743.1496 15624.397 3.1928399 636.43967 -226932.03 2843.7615 23.025673
      305 30.500001 -6000.8267 15605.599 3.1966859 609.29356 -226920.37 2722.4664 23.016435
      310 31.000001 -1063.0508 15589.326 3.2000228 587.11021 -226922.64 2623.346 23.008432
      315 31.500001 2193.2331 15576.216 3.2027161 569.89225 -226940.5 2546.4121 23.00198
      320 32.000001 2567.0546 15566.373 3.2047413 557.53042 -226975.38 2491.1765 22.997134
      325 32.500001 1335.9831 15559.193 3.20622 549.79307 -227027.8 2456.6042 22.993598
      330 33.000001 900.59352 15553.835 3.2073246 546.25734 -227096.9 2440.8057 22.990958
      335 33.500001 3261.0909 15549.861 3.2081443 546.28734 -227180.63 2440.9398 22.989
      340 34.000001 7394.1615 15547.588 3.2086133 549.29282 -227277.15 2454.369 22.98788
      345 34.500001 10393.395 15547.715 3.208587 554.89962 -227385.79 2479.4215 22.987943
      350 35.000001 10927.608 15550.55 3.2080021 562.94917 -227506.77 2515.3888 22.98934
      355 35.500001 10334.892 15555.745 3.2069309 573.24408 -227639.79 2561.3889 22.991899
      360 36.000001 10792.287 15562.757 3.2054858 585.4279 -227783.8 2615.829 22.995353
      365 36.500001 12554.247 15571.397 3.2037071 598.96032 -227937.06 2676.295 22.999608
etc

Arthur

Le 10/04/2015 21:06, Axel Kohlmeyer a �crit :

The paper I found is here: http://www.wag.caltech.edu/home/duin/papers/SiSiO2.pdf

If you scroll down, they report density as 4.29 g/cc. The other densities they report for quartz, cristobalite, and coesite I get really close to their values. That is why I know the FF is correct. But for stishovite, I have no idea how it remained steady at 4.29 g/cc.

I have changed the timestep, the frequency of updating the self-charges on the atoms (I usually just do it every time step, but tested other frequencies just to be sure), the number of atoms etc.

Also I have used both ReaxFF potentials by Kulkarni et al. and Fogarty et al. with similar results.

Ben

I am referring to Table 6 in the link I sent you. Here it is again: http://www.wag.caltech.edu/home/duin/papers/SiSiO2.pdf

Is there a track record of ReaxFF SiO being able to reproduce correctly
dynamics of stishovite? In the original paper (van Duin, J phys chem A,
2003), it looks like all that has been done for stishovite (and some
other polymorphs) are single point energy calculations at different
densities, to check the equation of state. I don't want to criticize
ReaxFF, which is a brilliant approach to reactive forcefields, but this
is a typical flaw often present in papers introducing new
parametrizations of forcefields. It is not enough to compute the energy
of a system, keeping its structural aspect constant and only varying its
density (isotropically in general). It's not enough either to only
perform energy minimization. The only correct way is to perform dynamics
at a specific pressure and temperature.

I second Arthur on the note about dynamics which is nothing but to say that
the approach/parametrization should follow more a force-matching method
like originally
proposed by Ercolessi.

I've tried to reproduce your results using reax/c, and using the FF
parameters file included in the supporting information of the J phys
chem C paper (Kulkarni, 2012). It's important to note that this is a
reparametrization of the original ReaxFF SiO forcefield, which is also a
possible explanation of why stishovite isn't correctly described.
Anyway, I'm getting something very similar to what you get:

Were the results of the paper produced using Lammps? (no access to
e-resources on this side right now) Van Duin has his own reaxff fortran
code which also uses an external file (like lammps) where some of the FF
params can be tuned. Sounds like Ben has contacted him already so perhaps
is a good idea to get the fortran code to run the same test before rolling
the dice that points to which side of the fence is where the mistake lives.
Carlos

I did indeed contact Dr. van Duin but am still awaiting a reply. I was never able to get a hold of the original 03 FF. Dr. van Duin said that Fogarty et al (2010). would be sufficient. And like I said, Kulkarni et al. (2013) produces similar results as well.

Ben

Stishovite is a high pressure, metastable silica phase. It is likely you cannot equilibrate it at zero pressure. You may also need smaller time step and much larger system sizes.

Ray

I disagree - yes, stishovite is metastable, but that does not imply that it cannot be equilibrated at ambient T/P. If it was the case, there would be no samples of stishovite found in nature at ambient thermodynamical conditions. As an other example we could not be able to model diamond carbon either, which is metastable too. Yes, it would be a good idea to test this on the original fortran version of ReaxFF. I believe that the results of the paper were produced using Van Duin’s fortran ReaxFF, if I’m not wrong ReaxFF wasn’t included in LAMMPS at this time (2003). Arthur

It might be the case that I cannot equilibrate at zero pressure for ReaxFF, but other potentials perform quite well, including the one you developed, COMB10.

Ben

I imagine you already tried plain energy minimization prior to the npt equilibration, right?

On another front, if you print out during your npt run the energy components of the reaxff FF (broken down into its individual terms) can you spot the issue being triggered by one particular term or all of them seem to contribute to the temperature increase? You may want to reduce the timestep a bit for this exercise.

Carlos