I am trying to simulate a Fe-Cr system (using EAM potential) with a single dumbbell interstitial defect and a system size of 2001 atoms at 1000K. I am tracking the defect using LAMMPS Voronoi occupation. I use the following (key) steps in my LAMMPS script:
Create a Fe lattice (system size is 2000 atoms) and minimize using ‘fix box/relax aniso 0’.
Introduce the dumbbell interstitial in place of an atom (now system size is 2001 atoms) and minimize the system again using ‘fix box/relax aniso 0’.
Introduce the Cr atoms (at.% 1) and minimize the system again using ‘fix box/relax aniso 0’.
Remove all the fixes and heat the system using ‘velocity all create 1000 56889 mom yes rot yes’ and apply the following fixes
fix LGV all langevin 1000 1000 1 56889
fix PRESS all press/berendsen aniso 0.0 0.0 500 couple none
fix NVE all nve
fix MOM all momentum 1 linear 1 1 1 rescale
run 100000
I use ‘fix MOM all momentum 1’ because I want to make sure the system does not drift at all during the heating since the reference system has already been created for tracking the defect when the lattice was created. And if the system moves/drifts, the defect cannot be tracked anymore.
I remove all the fixes and apply new ones as follows
fix NVE all nve
fix MOM all momentum 10000 linear 1 1 1 rescale
And run the system for 250 million timesteps and start tracking the defect. The system runs smoothly up until ~100 mill. timesteps, post which it drifts, and my results become inaccurate.
Is there a way to stop the drift throughout the simulation of 250 million timesteps?
You describe your simulation to humans instead of just sharing the input deck. Are you expecting anyone to run your simulation in their head based on how they interpret/remember LAMMPS commands? Well, good luck with the answers!
I can only make two educated guesses:
If a sample drifts, it can be due to poor conservation of energy due to incorrect settings (eg large timestep) or to excess kinetic energy due to the creation of defects.
Having a fix momentum in the production run is hardly a good idea. If the sample still drifts, increase the equilibration run by 1M steps.
Thanks for the honest response and the suggestions. I can try an equilibration run of 1M steps. The current timestep is 0.001 ps (default from ‘units metal’).
Thanks for sharing the guidelines, I will edit my first post accordingly.
About not sharing the input script, I figured it might be better to list the key steps with fixes (a pseudo code) as sharing the script seemed like, “Hey, tis mah problem, solve it! errr please.” Plus, the drift seems to appear after 100 million timesteps only, so for someone to recreate and observe the issue, they would need hours.
This is a key problem. It leaves no room for equilibration.
This creates a problem because it interferes with the dynamics and suppresses fluctuations.
Depending on the kind of system, it can be nearly impossible to suppress some (mild) drifting while the system is still equilibrating. Only after the system is in equilibrium you can assume that the system will not drift and any resulting drift is the result of accumulating numerical errors (mostly from using floating-point math).
While in your initial part you apply fix momentum too often, it could be applied more often in step 5, say every 1000 steps. The longer you wait, the larger the disruption when you remove and redistribute the energy from the center of mass momentum.
I also would experiment to replace fix langevin, fix press/berendsen and fix nve with fix npt.
One more option to consider would be using fix recenter. This will not remove any drift, but it can remove the motion of your region of interest. So it would go on top of everything.
If you are concerned about the centre of mass drifting then you need to use the zero yes option of fix langevin.
You may also want to consider removing the rescale keyword, since it introduces more floating point operations and thus reduces the accuracy of your intended momentum conservation.