Dear LAMMPS Users,
I am using the Core-Shell model to Simulate SrTiO3 using a potential that I found in literature (Speliarsky et al., Current Opinion in Solid State and Materials Science, 9, 107, 2005), this paper gives all the parameters for the Buckingham + Coulomb potential. When I run the simulation atoms start to “fly” outside the box while the temperature increases dramatically and then I get an error. This only happens at higher temperatures above 80K.
I tried debugging the script by:
-changing the timestep
-Different Pdamp and Tdamp
-Starting with different atomic coordinates (i also just took the restart file from previous runs)
-Different core shell mass ratio for oxygen
To me it looks like the oxygen atoms move too far and bond to some other core so in the simulations a core ends up with two shells (is this possible?) and thus the system explodes.
Is there something I am doing wrong or missing? I have also attached my LAMMPS input files in case someone wants to take a look.
Any advice would be much appreciated.
Best regards,
Ari
in.cooling (3.03 KB)
SrTiO3.lmp (1.77 KB)
I suggest you massively simplify your system to isolate the cause of the instability.
from your description, first item to look into would be the O-O core-shell bond.
to test this, you can set up a system with fixed size box with 1 O core and O shell atom and only the bond between them and only time integration on the shell.
that would allow you to validate the bond parameters from the literature (chances are you made a mistake when converting them) and you can find out what time step is required for a stable time integration of that bond alone. it may be smaller than what you currently have.
then i would add another atom and its shell and - again - apply the time integration to the shell only and adjust the timestep if needed (probably not, considering that that O atoms are giving you trouble).
then add the last type of atom and its shell and repeat.
now set up a periodic system with the pair style only and no shells and determine whether this can run without crashes and with good energy conservation. i would start this with fixed volume and only switch to isotropic fix npt after it has releaxed, then anisotropic and tri only if really needed. allowing the system to tilt can lead to problems, since a maximally tilted box formally has a very low energy but is unphysical and unstable. after tweaking the settings for this setup to get good energy conservation, you have now validated parameters for the non-bonded and bonded terms and also optimal time steps for each.
the final step would be to put all together, you possibly want to use r-RESPA timestepping so you can use a (much?) smaller timestep for the core-shell bonds only. with a protocol like this, it should be easier to narrow down where the problem is and then investigate closer how to resolve it.
axel.