reax and reax/c minimization question

Hi all,

I ran minimizations using the reax/rdx examples and found them being very different. While the energy minimizes to similar value, the forces are by more than an order of magnitude. The log files are attached.

The pair_style reax minimizes to force ~ 24 (real units) and reaxc to 0.5. It seems the reaxc has a more reasonable results given the purpose of minimization.

The same phenomena happend to my other systems. While reax/c minimizes the force better, it crashes in production run when I tried to increase the timestep. Reax does not crash with larger timestep. So I was really hoping I could use the two interchangeably.

Kan-Ju

log.rdx_reax (1.95 KB)

log.rdx_reaxc (2.83 KB)

Kan-Ju,

For some reason that is also not clear to me, minimizer cg does not seem to work well with pair_style reax for this test case. From the manual, hftn is suggested as an alternative when cg seems to perform poorly.

Using hftn with pair_style reax, I obtained the following, which does not look bad.

Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-1884.31235962 -2039.60503845 -2039.61631199
Force two-norm initial, final = 1039.11 0.0487185
Force max component initial, final = 346.372 0.0266989
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 20 909

Ray

Hi Ray,

Thanks for the response. You are right, changing to hftn does help for the rdx reax example. I did the tests for my systems (for both reax and reax/c), but didn’t find this improvement.

I’ve done hftn and cg with backtrack/quadratic/forcezero on both reax and reax/c for alpha quartz (72 atoms). The log files and the lammps inputs are attached.

With different min_style and min_modify, the force varies 4 to 1575. The smallest force I got out is 5.0 (reaxc/cg_quadratic) and 17.6 (reax/cg_forcezero).

I’d like to know if I could mix the minimized energies obtained from reax and reax/c. For example, to obtain the formation energy from SiO2 and O2 minimized from reax/c, and Si minimized from reax.

Thanks.

Kan-Ju

logfiles.zip (10.6 KB)

SiO2_minimization.zip (13.1 KB)

Kan-Ju,

Sorry, but I will not be able to find time to go through all your files. My suggestion is to identify the one pair_style and minimizer that gives you an overall acceptable and reasonable convergence and stick with it. It is usually advised against mix-and-match.

You can not rule out the possibility of falling into local minima using minimizations, and if the purpose is to calculate heats of formation you may end up with bogus results. You may try using a low temperature annealing using NPT, and block-average the equilibrium portion of the potential energy evolution.

Ray

Hi,

I totally understand. Thanks for helping.

The ideal case would be one setting for everything. From the systems I’ve tested (from primitive to 10x10x10), the performance varies sometimes by 2 orders of magnitude. My thought is that since all the minimization methods are ways to minimize the energy and the final force/energy are evaluated using the same potential, it might be alright to use any approach to reach the optimal status. I have little knowledge on those algorithms so perhaps I am wrong. Despite the huge force difference, the energies come out to be the same.

One thing that confuses me is the outcome of reax and reax/c. I used the minimized structure from reax (E=-2363), and perform single point calculation (run 0) using reax/c. The reax/c gives a much higher energy (-2239), while its minimized energy is at -2349. Is this difference meaning the potentials are evaluated differently in the two packages?

I was using reax/c for all my systems until it keeps crashing when I try to increase the timestep. I’ve let the system to run at 0.04 fs for half million timesteps, but couldn’t move up to 0.05 fs (the goal is 0.25). I started to test reax after I found that reax runs successfully for 0.05 fs.

The error message for reax/c for timestep 0.05:

Hi Kan-Ju,

Pair_style reax/c is typically 2 or more times faster than pair_style reax. You can take a look at the benchmark page: http://lammps.sandia.gov/bench.html#potentials. Different minimization methods use different algorithm to find the minima - some gradient based and some damped dynamics, so it is expected that each minimizer gives slightly different result.

The stand-alone ReaxFF code from Cal Tech has a control file that controls a wide variety of variables/parameters that are not in the ffield.reax files. These control settings are hard-coded in reax, but are free for users to change in reax/c (via the control file). It is unlikely that the results differ because of these settings, but it is one possibility. Were the same minimizer used for both runs? Generally reax and reax/c give very comparable results.

The error occurred because the number of non-bonded interactions increased dramatically during the run, and it ran out its allocated array size. You can use the safezone and mincap keywords to increase the buffer, which should help in this case.

Ray

Hi Ray,

Thanks. Increasing the safezone to 1.6 and mincap to 100 allows the timestep 0.05 to run without any problem.