Dear all,
I am trying to apply a potential that I found in the literature (T.S. Bush et al, 1994) for Ga2O3. The paper gives the parameters for a Coulomb-Buckingham potential and a core-shell model with harmonic bonds. Below is my attempt to use this model for a beta-Ga2O3 crystal, starting from the adiabatic core shell example. When I run the script and visualize the xyz file in Ovito, what I see is that the shells are being “ejected” from the cores (see picture) as the temperature increases: it reaches over 2000000K at the last step! The temperature also increases in the example file, but the cores and shells seem to remain bond, and the temperature eventually decreases after reaching a peak at about 500K.
I tried debugging the script by myself, but I am not sure how to interpret this effect. However here is what I observed:
-
when I increase dt, the effect that I describe above happens faster and stronger, and bonds are quickly lost. I suspect that this happens when a shell or core crosses a boundary and finds itself on the other side of the simulation box, although I am not quite sure about that.
-
I don’t have parameters for the Ga-Ga interaction, as opposed to the Na-Cl example. Moreover the parameters given in the paper I mentioned above are an order of magnitude lower to the ones used in the Na-Cl example. Maybe the harmonic bonds are not strong enough to keep the cores and shells together?
-
the line with special_bonds, which I intended for coulombically shielding the shells from the cores, makes no difference at all in the result, as if the shells were already shielded before: why is that?
Any idea what I am doing wrong? Thank you very much in advance for your answers and ideas.
Julien Barrat.
------------------------ Initialization -------------------------------
units metal
atom_style full
boundary p p p
------------------------ Atom Definition -------------------------------
fix csinfo all property/atom i_CSID
read_data Ga2O3.data fix csinfo NULL CS-Info # this file contains a “rectangularized” unit cell of beta-Ga2O3
unfix csinfo
replicate 4 4 4
group cores type 1 2
group shells type 3 4
neighbor 2.0 bin
comm_modify vel yes
------------------------ Force Fields -------------------------------
kspace_style ewald 1.0e-6
pair_style buck/coul/long/cs 12.0 # cutoff
pair_coeff * * 0.0 1.0 0.0
pair_coeff 3 4 2339.776 0.2742 0.0 # Ga-O (shells): A, rho, C
pair_coeff 4 4 25.41 0.6937 32.32 # O-O (shells): A, rho, C
bond_style harmonic
bond_coeff 1 10.265 0.0 # k=20.53 eV
bond_coeff 2 10.265 0.0 # k=20.53 eV
special_bonds coul 0.0 0.0 0.0
------------------------ Equilibration Run -------------------------------
reset_timestep 0
thermo 100
thermo_style custom step etotal pe ke temp press
compute CStemp all temp/cs cores shells
compute thermo_press_lmp all pressure thermo_temp # press for correct kinetic scalar
velocity bias option
velocity all create 300 134 dist gaussian mom yes rot no bias yes temp CStemp
velocity all scale 300 temp CStemp
thermostating using the core/shell decoupling
fix thermoberendsen all temp/berendsen 300 300 0.4
fix nve all nve
fix_modify thermoberendsen temp CStemp
dump id all xyz 10 potential.xyz
timestep 0.000001
run 10000
