Dear all,
I ran NVT ensemble on a box of Ga2O3 at 3000 K using buckingham fordefield. Initially the temperature was correct, which is around 3000 K, then T became zero. I am wondering what is the reason and how to solve it. I posted my input file below:
variable name string ga2o3 # Molecule id
variable T equal 3000 # Simulated Temperature [K]
variable tempname equal 3000 # Temperature Name in Files
variable dt equal 0.0005 # timestep
variable runnumber string 02 # Output File Name
variable prunnumber string 01 # Input File Name (manually change initial in first run)
units real
dimension 3
boundary p p p
atom_style charge
read_restart ga2o3-2.lammps-run01
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes page 100000
kspace_style pppm 1.0e-5
pair_style buck/coul/long 10.0 12.0
timestep ${dt}
fix 1 all nvt temp 3000.0 3000.0 0.1
thermo_style custom step temp pe etotal press density vol
thermo 1000
run 1000000
i am wondering where is wrong. Thank you!
Have you visualized your trajectory?
Are there any warnings or other unusual output in your log file?
Not really. The following is the output which showed the abrupt change in temperature during the simulation:
Step Temp PotEng KinEng TotEng Press Volume
9700 3239.6951 -21233012 250297.6 -20982715 -544787.53 271517.48
9800 61.667413 -21240715 4764.4006 -21235950 -587389.93 271517.48
9900 189.26789 -21262338 14622.764 -21247716 -588567.64 271517.48
10000 14755.165 -21307438 1139978.4 -20167459 -404373.56 271517.48
10100 1470.4416 -21339779 113605.76 -21226173 -580465.25 271517.48
10200 0.00026033243 -21261243 0.020113183 -21261243 -591780.54 271517.48
10300 0.00026581614 -21261464 0.020536853 -21261464 -591807.53 271517.48
10400 0.00027146286 -21261687 0.020973116 -21261687 -591834.74 271517.48
10500 0.00027727937 -21261913 0.021422497 -21261913 -591862.2 271517.48
10600 0.00028327282 -21262141 0.021885549 -21262141 -591889.9 271517.48
10700 0.00028945076 -21262372 0.022362854 -21262372 -591917.84 271517.48
10800 0.00029582116 -21262606 0.022855029 -21262606 -591946.05 271517.48
10900 0.00030239246 -21262842 0.023362725 -21262842 -591974.52 271517.48
11000 0.00030917358 -21263082 0.023886632 -21263082 -592003.27 271517.48
11100 0.00031617395 -21263324 0.024427478 -21263324 -592032.29 271517.48
11200 0.0003234036 -21263569 0.024986038 -21263569 -592061.6 271517.48
11300 0.00033087312 -21263818 0.02556313 -21263818 -592091.22 271517.48
11400 0.00033859377 -21264069 0.026159625 -21264069 -592121.13 271517.48
11500 0.00034657753 -21264324 0.026776447 -21264324 -592151.37 271517.48
11600 0.00035483711 -21264583 0.027414579 -21264583 -592181.92 271517.48
11700 0.00036338605 -21264845 0.028075066 -21264845 -592212.81 271517.48
11800 0.00037223874 -21265111 0.028759022 -21265111 -592244.05 271517.48
11900 0.00038141058 -21265380 0.029467635 -21265380 -592275.65 271517.48
Just look at these numbers. Between steps 9700 and 10000 something very drastic has to have happened to your system. It is a well known property of the Buckingham force field formulation that they are prone to a “Coulomb collapse” due to lack of repulsion at short range. Please also note that your pressure is very unusual for a high temperature condensed system. This suggests that either there is a problem with the force field parameters, or geometry, or both. A careful inspection and visualization of the trajectory is the first step to diagnose this.
Watching the trj file using vmd, I can see that before 9700 steps, the crystal was stable. At that step, I saw one stretching in the box. So one atom did an abrupt movement.
The system has 25920 atoms. The dt = 0.001 in my test, which is already very small. I also did one test by reducing it to a half, but after double the number of steps, I still saw the temperature decreasing to zero. Besides that, I also did the test reducing kspace pppm to be 10e-6, instead of pppm to be 10e-5. After running the simulations for sometime, I can see the temperature dropped and approached to zero. I am wondering what else i can try at this moment to solve the issue.
With such a small timestep you are not sampling much of the configurational space. My guess is that the initial structure is out of equilibrium and that by firing your sample at 3000K, one or more atoms overcome the Buckingham energy barrier and the simulation fails. Or better it doesn’t, because you are using such a small timestep. I would do the following:
- Print a
thermo output which is more informative, e.g.: thermo_style custom step etotal evdwl evdwl ecoul elong pe ke temp press
That will help you get a better physical picture.
- Start by equilibrating the sample at T = 300 K, dt = 1 fs, and the recommended coupling constant for the thermostat, e.g.
fix 1 all nvt temp $T $T $(dt*100)
(also, why don’t you use the variable $T if you define it?)
- Then switch to NPT and increase the temperature from 300 to 3000 over a few nanoseconds (stick to dt = 1 fs).
Probably the simulation will fail at some point, due to the Buckingham potential not being suitable for high-temperature melts. A trick I used to anneal fused silica was to fit the Buckingham potential with a Lennard-Jones potential, which has a hard repulsive barrier, and use the modified force field only at high temperature. When you cool the sample below a safe temperature (which you will learn from point 3), you can switch to the original force field and equilibrate well before extracting any meaningful information from the simulation.
Just my 2 cents.
1 Like