ReaxFF Efficiency

Greetings LAMMPSters!

Recently I am trying to get some production data in LAMMPS and it seems to be quite slow. I am seeking advice on how to speed up my calculations. The reactions I am simulating have a fairly high activation energy (~10-15 kcal/mol) and take place in solvent, so it is imperative to have the fastest computation times possible.

Are there any tips or tricks to speed up your computations using ReaxFF?

Is there anything obviously wrong or anything that would considerably slow down computations in the following script:

REAXFF Simulation

units real
atom_style charge
atom_modify map array sort 1000 0.0
read_data data.system
read_dump dump.system 401000 x y z box yes

pair_style reax/c lmp_control
pair_coeff * * ffield.reax.chosi C H O Hc C O Si Hc H O
neighbor 2 bin
neigh_modify every 10 delay 0 check yes
fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c
variable Density equal mass(all)*1.66053892/vol

thermo_style custom cpu step time etotal pe temp press v_Density
thermo 500
dump 1 all atom 500 dump.restart

fix 3 all restrain bond 4609 4610 50.0 50.0 1.52
fix_modify 3 energy yes
fix 4 all restrain bond 4611 4620 50.0 50.0 1.42
fix_modify 4 energy yes
fix 5 all restrain bond 4612 4617 50.0 50.0 1.42
fix_modify 5 energy yes
fix 6 all restrain bond 4613 4614 50.0 50.0 1.42
fix_modify 6 energy yes

Lots of bond restraints

fix 4992 all restrain bond 4605 4606 50.0 50.0 1.10
fix_modify 4992 energy yes
fix 4993 all restrain bond 4605 4607 50.0 50.0 1.09
fix_modify 4993 energy yes
fix 4994 all restrain bond 4605 4608 50.0 50.0 1.09
fix_modify 4994 energy yes

fix 4995 all nvt temp 2000.0 2000.0 100.0
fix 4996 all momentum 1000 linear 1 1 1
fix 4998 all reax/c/bonds 1000 bonds.reaxc

timestep 0.25
run 20000000

Currently, I am simulating at high temperatures (1500K-2500K) to even see these reactions happen, and I have been using fix restrain in order to prevent decomposition reactions at high temperatures.

I have tried editing the following molecular dynamics parameters and obtain a small or negligible speed benefit:
-fix qeq/reax convergence cutoff
-fix qeq/reax Nevery
-neighborlist Nevery
-In lmp_control, granularity of long range tabulation (currently using 1000 points)

I have also tried benchmarking my version of lammps and find that I am about 20-30% slower, which is expected as I am using 6 y/o hardware.

On scaling my system I find shorter cpu-time per timestep.
With a system of 8748 atoms. The cpu time for 100 steps:
4 processors: 45s

8 processors: 27s
16 processors: 17s
32 processors: 11s
64 processors: 8.3s

I have also tried using different compilers and optimization levels with lammps, but I only found 3-5% difference in each case.

Best,
Josh Deetz
PhD Candidate Student
Chemical Engineering
University of California, Davis

Comments below, thanks.

Ray

Greetings LAMMPSters!

Recently I am trying to get some production data in LAMMPS and it seems to
be quite slow. I am seeking advice on how to speed up my calculations. The

Slow and fast are all relative. From your scaling data at the end, it
does not seem slow at all that you can run 100 steps in 8.3 seconds
with 64 MPI processes (that is equivalent to ~260 ps/day with a 0.25
fs timestep). In fact, this can be considered fast for ReaxFF.

Comments below, thanks.

Ray

Greetings LAMMPSters!

Recently I am trying to get some production data in LAMMPS and it seems to
be quite slow. I am seeking advice on how to speed up my calculations. The

Slow and fast are all relative. From your scaling data at the end, it
does not seem slow at all that you can run 100 steps in 8.3 seconds
with 64 MPI processes (that is equivalent to ~260 ps/day with a 0.25
fs timestep). In fact, this can be considered fast for ReaxFF.

reactions I am simulating have a fairly high activation energy (~10-15
kcal/mol) and take place in solvent, so it is imperative to have the fastest
computation times possible.

Are there any tips or tricks to speed up your computations using ReaxFF?

If you think hydrogen bonds are not important in your system, you can
decrease its cutoff using control settings - this usually speeds
calculations up. You can also decrease fix qeq/reax cutoff.

Is there anything obviously wrong or anything that would considerably slow
down computations in the following script:

# REAXFF Simulation
units real
atom_style charge
atom_modify map array sort 1000 0.0

Sorting would slow things down.

read_data data.system
read_dump dump.system 401000 x y z box yes

pair_style reax/c lmp_control
pair_coeff * * ffield.reax.chosi C H O Hc C O Si Hc H O

Does ReaxFF have an element called Hc?

neighbor 2 bin
neigh_modify every 10 delay 0 check yes
fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c

You can decrease the "10.0", or the cutoff to possibly 7 or 8 to speed
things up a little bit.

variable Density equal mass(all)*1.66053892/vol
thermo_style custom cpu step time etotal pe temp press v_Density
thermo 500
dump 1 all atom 500 dump.restart

fix 3 all restrain bond 4609 4610 50.0 50.0 1.52
fix_modify 3 energy yes
fix 4 all restrain bond 4611 4620 50.0 50.0 1.42
fix_modify 4 energy yes
fix 5 all restrain bond 4612 4617 50.0 50.0 1.42
fix_modify 5 energy yes
fix 6 all restrain bond 4613 4614 50.0 50.0 1.42
fix_modify 6 energy yes
...
Lots of bond restraints
...
fix 4992 all restrain bond 4605 4606 50.0 50.0 1.10
fix_modify 4992 energy yes
fix 4993 all restrain bond 4605 4607 50.0 50.0 1.09
fix_modify 4993 energy yes
fix 4994 all restrain bond 4605 4608 50.0 50.0 1.09
fix_modify 4994 energy yes

fix 4995 all nvt temp 2000.0 2000.0 100.0
fix 4996 all momentum 1000 linear 1 1 1
fix 4998 all reax/c/bonds 1000 bonds.reaxc

timestep 0.25
run 20000000

Currently, I am simulating at high temperatures (1500K-2500K) to even see
these reactions happen, and I have been using fix restrain in order to
prevent decomposition reactions at high temperatures.

These fix restrain bond commands are against the ReaxFF spirit. If
something other than what you want to see decomposes at high
temperatures, then this parameterization is not suitable for what you
want to model.

I have tried editing the following molecular dynamics parameters and obtain
a small or negligible speed benefit:
-fix qeq/reax convergence cutoff

There is no such thing as "convergence cutoff". Either the Taper
function cutoff or the tolerance.

-fix qeq/reax Nevery

Changing this frequency is not recommended.

Ray, thank you for your recommendations and taking the time to read my long message!

Hc is for methyl hydrogens, H is for hydroxyls and water. I did this because the vdW parameters should be different for Hc and H.

I have one other question: will specifying atom types with “pair_coeff * * ffield.reax.chosi C H O Hc C O Si Hc H O” influence the speed of the simulation? Do you think condensing this to “pair_coeff * * ffield.reax.chosi C H O Si Hc” would be any faster?

Thanks,

Josh

Ray, thank you for your recommendations and taking the time to read my long
message!

Hc is for methyl hydrogens, H is for hydroxyls and water. I did this because
the vdW parameters should be different for Hc and H.

I have one other question: will specifying atom types with "pair_coeff
* * ffield.reax.chosi C H O Hc C O Si Hc H O" influence the speed of the
simulation? Do you think condensing this to "pair_coeff * *
ffield.reax.chosi C H O Si Hc" would be any faster?

You are welcome. This should not affect the efficiency, but you can
run a quick test to verify that.

Ray

Hi Joshua

Are you simulating zeolite synthesis or catalytic reactions? I am very interested because we are working on such things.

As Ray said, ReaxFF is intrinsically slow compared to classical force field. There is few acceleration techniques for ReaxFF as far as I know, except GPU, which significantly reduces the cost of QEq. See http://dx.doi.org/10.1016/j.jmgm.2013.02.001

Simulation at high temperature is simple and effective, but also raises problems as the distribution is altered.

Our previous work use Simulated Tempering. http://dx.doi.org/10.1021/jp300221j

And now we are focusing on other enhanced sampling methods such as REMD and accelerated MD. For them to be most efficient, the implementation is not straightforward.

You may try Parallel Replica Dynamics, which has recently been validated by Adri, although it is designed for small systems. http://dx.doi.org/10.1021/jz4019223

Your use of two Hydrogen types is especially interesting. Could you explain the reason for me? And do you allow methyl H and hydroxyl H to exchange?

Regards,

Francis

Greetings LAMMPSters!

Recently I am trying to get some production data in LAMMPS and it seems to
be quite slow. I am seeking advice on how to speed up my calculations. The
reactions I am simulating have a fairly high activation energy (~10-15
kcal/mol) and take place in solvent, so it is imperative to have the fastest
computation times possible.

sorry, but i disagree on this statement. you don't need to run faster,
but you need to run smarter. why are you using such a brute force
approach to get the system what you want to do? what is the scientific
value going to be? you don't say what kind of information it is that
you are after, but it seems you could just as well, make up a fake
trajectory that shows the reaction you want to see, which means you
are making a movie, but not gathering scientific evidence.

there is a whole zoo of methods around that provide enhanced sampling
in order to better deal with rare events and LAMMPS supports several
of them. why not define your desired reaction pathway(s) as collective
variables and apply a suitable biasing method?

axel.

Axel,

Thank you for your reply. In fact thank you for responding to all of my messages in the past, I owe a great deal to you and others on the mailing list.

I will provide you with more details about my project. My system has three molecules types in the initial configuration (A=200 molecules, B=1600 molecules, C=1600 molecules). All three of these molecules can react with eachother. My initial plan was to use parallel tempering with this system using the replica package in order to drive the system down to room temperature. Is there any other method for accelerating such a system?

-Josh

Axel,
Thank you for your reply. In fact thank you for responding to all of my
messages in the past, I owe a great deal to you and others on the mailing
list.

I will provide you with more details about my project. My system has three
molecules types in the initial configuration (A=200 molecules, B=1600
molecules, C=1600 molecules). All three of these molecules can react with
eachother. My initial plan was to use parallel tempering with this system
using the replica package in order to drive the system down to room
temperature. Is there any other method for accelerating such a system?

accelerate to get what information out of the simulation?

please keep in mind that due to the restrictions of time and length
scales it is almost *never* a solution to just stick some molecules
into a box and then run a simulation and then look at the outcome,
particularly if you have a complex system with many alternatives.

my guess(!) would be that MD may not even be the right method for what
you want and running some kinetic monte carlo simulations could be a
more promising path to follow.

axel.