Question on Suitable timesteps for reax vs eam forcefield

I am running a tension test on a set of atoms and trying different forcefields to test their suitability with producing an accurate stress strain curve. I began with the lammps tutorial simulation which uses the EAM forcefield and a timestep of 0.001 ps.

Next, I did this exact same script but did it with a reaxff forcefield. While I made the units “real” I forgot to switch the timestep, so it was 0.001 fs. However, the modulus was accurate. Ideally I want to compare simulations of the same length, meaning decrease strain rate and run longer. However, what is a recommended timestep for reaxff? Is 0.001fs too small? What about strain rate?

Thanks!

Liam

I am running a tension test on a set of atoms and trying different forcefields to test their suitability with producing an accurate stress strain curve. I began with the lammps tutorial simulation which uses the EAM forcefield and a timestep of 0.001 ps.

Next, I did this exact same script but did it with a reaxff forcefield. While I made the units “real” I forgot to switch the timestep, so it was 0.001 fs. However, the modulus was accurate. Ideally I want to compare simulations of the same length, meaning decrease strain rate and run longer. However, what is a recommended timestep for reaxff? Is 0.001fs too small? What about strain rate?

​the optimal time step depends on several properties, for example temperature, particle mass, steepness of the potential, density.
if you run the same system with different potentials, those should run at mostly the same settings. however, in reaxff you ​have charge equilibration, and that usually requires a somewhat shorter time step. what is optimal, you can figure out empirically. 1fs may be too aggressive but could perhaps work well enough for your particular system and 0.1fs would be a conservative choice. the way to detect this is to start from a longer time step and reduce it until results you are interested in don’t change significantly anymore. or do it the other way around.

axel.

Thanks for the quick response, I am going to experiment a bit the weekend. Where I assume I must also be careful is in the strain rate. If I switch the time step I will need to also need to be careful with strain rate. I am thinking I am best keeping strain rate the same for each simulation (ie changing it to make the units match). Right now for my EAM the strain rate was .01/ps (as per tutorial) but my reax is a much shorter simulation with a strain rate of 0.01/FS. Interestingly though the modulus for reax and EAM when set up like this both make sense.

Liam

I run a lot of solid state carbon simulations with ReaxFF and typically use time steps between 0.05-0.2 fs. Within the ReaxFF framework, the various parameter sets can require substantially different timesteps depending on the curvature of the energy profile. For example, with carbon, the CHO parameters of Chenoweth et al. require a shorter timestep than more recent C-2013 parameters of Srinivasan et al. due to a sharp change in the energy profile around the bond equilibrium distance in the CHO parameters. I’ve found low strains to be less sensitive to timestep than larger strains. Specifically, I’ve seen larger timesteps cause premature fracture, with no other obvious problems with the simulations.

-Ben

Thanks for the advice.

I have realized that when switching to reaxff I obviously must also ensure my strain rate stays consistent. However, I am coming across a very odd error on my script. When I change the timestep to 0.25 the lammps executable crashes once the equilibration starts. If I make it 0.025 it works, 1 it works, but if it is zero followed by a decimal and then any number it crashes? What in my script could be causing this or is there an error I am not aware of. 0.25fs seems standard for reaxff so I am confused as to why it isn’t working. Is it a windows issue?

ironscriptunitfix.txt (3.16 KB)

Thanks for the advice.

I have realized that when switching to reaxff I obviously must also ensure my strain rate stays consistent. However, I am coming across a very odd error on my script. When I change the timestep to 0.25 the lammps executable crashes once the equilibration starts. If I make it 0.025 it works, 1 it works, but if it is zero followed by a decimal and then any number it crashes? What in my script could be causing this or is there an error I am not aware of. 0.25fs seems standard for reaxff so I am confused as to why it isn’t working. Is it a windows issue?

you have a ​bad choice of Pdamp.​ it is far too small, which results in the system reacting too quickly to the total pressure with changing the volume and then creating a far too dense system, which will lead to all kinds of problems. you should not use the drag keyword unless your system shows unwanted oscillations. try a Pdamp of 10ps or 100.0ps.

as explained in the fix npt documentation, there is a relation between the size of the time step and the choice of Tdamp and Pdamp.

axel.

I must have missed the boxes in the npt documentation about this. Following the manual it says for T damp to be 100 time steps, so 25 fs. For P damp it says 1000 time steps so 250, however I also tried a Pdamp of 10 ps as you wrote above, or 10000fs. In all cases it came with an error saying fix qeq convergence failed at step 2 or 3. I also tried some other Pdamp, Tdamp combinations and got an error of not enough space for bonds.

Having done some reading about this it says it may be something due to my initial geometry, but what I am making is extremely basic. Simple unit cell. Again script is attached.

What I find most interesting though is the fact that if I run my initial script: 0.001 timestep, run for 2000 with a strain rate 0.01/fs, ran for a total of only two fs, produced a perfect stress strain curve. Maybe this is some lucky garbage in garbage out, but it only produced an accurate model when I added a fixqeq etc. In this script Tdamp and Pdamp were both 1 fs.

I am running something this basic so I can learn more about the program and how the forcefields work before I get into a more complex structure. Thanks for helping me out.

ironscriptunitfix.txt (3.18 KB)

you have exhausted my good will. this is all off-topic for the list, since this has nothing to do with LAMMPS, but is all about you learning MD, and doing this with a rather complex and difficult to use force field. you (or your adviser/supervisor) will have to find a suitable tutor for you. this mailing list is not the place.

axel.

Ok well I appreciate your help. My confusion is really just due to the fact that after following all the proper protocols I am now getting errors I have never seen before. Yet when I knew less and had the small time step I got “expected results”

Also if you are able to help me with this last issue I think it would allow me to run on my own from there. This site has been very helpful and I feel much more comfortable with the other forcefields. These errors are just new and occurring after I modify based on the manual (which is of course the best source).

Also if you are able to help me with this last issue I think it would allow me to run on my own from there.

​i don’t think so. you’ve made many errors so far with order of magnitudes, also you don’t seem to be understanding the time scales of ongoing and overlapping processes well enough, and finally, you seem to be too much hung up on having a single cause and a single explanation being the reason for something working or not working. that is not the case in MD simulations, especially with complex force fields like ReaxFF (you seem to be significantly underestimating its complexity).

besides, reviewing and advising you on your work is the job of your adviser. i don’t have the time to check your as thoroughly as needed. so far, there were always significant and easy to spot blunders. going beyond that requires much more time, care and effort.

like i was hinting at before, i am willing to give people that are new a certain amount of effort beyond what is on-topic for this mailing list to get them pointed into the right direction, but that is finite and you have exhausted this amount.

This site has been very helpful and I feel much more comfortable with the other forcefields. These errors are just new and occurring after I modify based on the manual (which is of course the best source).

​you just got lucky so far. mostly because you stayed far away from problematic areas and spent (wasted?) a lot of cpu time instead.​

​axel.​

Ok I understand and thank you for all your contributions. Perhaps the fact my modulus was accurate despite the low timestep was a bit lucky.

I will dig deeper to see why I am getting these errors. After searching the message boards I am seeing that perhaps due to my small system there is too many ghost atoms. I think my tdamp and pdamp are now ok, so I am thinking it is due to the fix qeq. I am basically learning this myself with the goal of mechanical properties from a stress strain curve.

Time learning it not time wasted! So I will keep at it :).

Cheers!

Just as a follow up for anyone following the thread, my tdamp was too low for y and z. When they were changed it ran and produced a good agreement!

If anyone has an opinion as to why when I ran with a timestep of 0.001 fs with t damp and p damp of 1 fs it worked I would really appreciate the thought. Even with this low timestep it ran and produced a very good elastic modulus.

Liam

Just as a follow up for anyone following the thread, my tdamp was too low for y and z. When they were changed it ran and produced a good agreement!

If anyone has an opinion as to why when I ran with a timestep of 0.001 fs with t damp and p damp of 1 fs it worked I would really appreciate the thought. Even with this low timestep it ran and produced a very good elastic modulus.

​as the documentation explains, tdamp and pdamp should be chosen differently​ according to the time step. if you have a 0.001 fs time step, 1fs is 1000x the time step and that is usually acceptable. if you have a time step of 0.1fs or 1fs, having 1fs for tdamp or pdamp is too small. since the pressure fluctuates much more for typical condensed systems than temperature, pdamp usually should be chosen larger than tdamp, at least 10x larger. the factors given in the documentation are for typical condensed systems with low compressibility. if this changes (e.g. very high temperature or very low density), also the optimal choice of tdamp and pdamp changes.

BTW: for figuring out the optimal choice of time step, you need to make test runs with fix nve and no thermostat/barostat (after equilibration) to see how well total energy is conserved.

to understand, why different choices for time step require different choices for pdamp/tdamp you have to understand nose-hoover chain thermostats, where you couple the system property (acceleration, either of atoms or the cell change) to chains of couple harmonic oscillators. you have to avoid strong coupling and resonances. tdamp and pdamp implicitly determine the characteristic frequency of the nose-hoover chains in combination with the time step as such.

in addition to this, problems always arise more easily when you get close to the maximum size of a time step where time integration is still stable.

​axel.​

Grab a textbook on real stress-strain curves of materials. If you ran with 0.001 fs and 2000 steps, you have only run 2 fs – that is not going to strain your system by a meaningful amount. Therefore you were sampling the very beginning portion of the elastic regime – hence the “perfect” stress-strain curve. The results obtained are therefore quite meaningless.

Ok that makes sense, thanks for the response. You are right I am certainly no where near the plastic regime. However, based on my strain rate I did strain 2 percent. The 0.001fs script strained at 0.01 per fs, meaning 2% after 2000. Based on my knowledge of solids, strains greater than 10% violate assumptions and staying below 5 % is advisable. So I while I agree simply showing the elastic modulus is what is should be is meaningless, am I correct in saying it does show that the reaxff forcefield is able to predict the mechanical property elastic modulus?

My surprise was that despite the high strain rate (0.01/fs) I still got a predicted result.

As a follow up to Rays comments, I took the very small timestep script (0.001 fs) run 5000 iterations to a total strain of 2% and did an nve energy check. After equilibrating for 10000 steps to 300k, the KE + PE was consistent. Moreover, when this script was run for stress strain it predicting a modulus of 203 (in linear portion of the stress strain curve). So perhaps the smaller timestep is also suitable.

Thanks all for your help, I have learnt a lot in this exercise.