Compute cross-interaction LJ parameters between materials and non-reactive water in a reaxff involved hybrid pair_style

Hi all,

I would like to use the OPLS force field to describe the bulk materials and reaxff to describe water and hydroxide ions. So I use the hybrid pair_style and explicitly specify the LJ parameters, similar to what this post did https://matsci.org/t/proper-usage-of-hybrid-usage-for-reaxff-and-gaff-force-field/34025. My question is: it is physically acceptable to use the parameter of Ow and Hw from non-reactive water model (spc/tip3p/tip4p) to compute the cross-interaction parameters between bulk materials and water based on mixing rule assuming non-reactive interaction between them?

Since my system is large (above 20000 atoms), I currently can not afford to model all atom types with reaxff and also there is no good parameterisation for the system. My primary objective is to observe how hydroxide interacts with water where hydrogen bonding is involved.

I wonder if the method mentioned above is a good alternative to deal with my system. Any suggestions are highly appreciated!

Best,

Eiji

It may be easier to answer your question if you can tell us more about your system. Depending on whether you’re modelling water around a protein, water mixed with ethanol, or water on the surface of an electrode, the simulation techniques available to you might be very different.

However, mixing ReaxFF with other force fields is almost always fundamentally a bad idea. ReaxFF is parametrised for short-range electrostatics with charge equilibration; most other MD force-fields (including OPLS) use fixed charges, and must thus include polarisation effects mostly in the LJ parameters. These two approaches are basically different and I would expect a “hybrid” simulation to be wrong unless proven otherwise.

Hi srtee,

Thank you for your reply. I am not sure how detailed I need to describe my system. Anyway, the system contained polymers tethered with ion-exchange functional group (sulfonate group), water, and KOH electrolyte, where OPLS-AA is used for bulk materials, tip4p is used for water, parameters for K+ and OH- are taken from literatures. My main goal is to see how hydroxide interacts with water to form hydrogen-bonding network within this polymer matrix. This is why I plan to use reaxff to account for the bond breaking and formation between OH and water. Since the system is a little bit large (may reach up to 30000 atoms), it is too expensive to model all atoms with reaxff. So I am looking for a feasible way to model the hydroxide transport.

Currently one way I can think of is to reduce the size of the sytem, maybe to a few thousand atoms. So the modelling and simualtion can be more convenient even using reaxff or ab-initio MD techinques. But I still wonder if there is any possible method to achieve this without changing too much on the original system (a relatively large system).

Best,

Eiji

Eiji,

a couple of suggestions:

  • if you don’t have enough computational power, you can always go an apply for CPU time on a national (or international?) supercomputing center. You may need a helping hand from your supervisor with the paperwork. What will work in your favor is that the ReaxFF pair style has been ported to KOKKOS and thus can be used effectively with GPUs. There is far more competition for CPU resources and the big centers are always looking for ambitious projects that require large amounts of time and can make good use of GPUs. That may allow running with ReaxFF throught after all. But I would be concerned of the usefulness of using ReaxFF for your system at all. People have a tendency to overestimate the applicability and accuracy of ReaxFF for systems where electrostatic and the balance between coulomb and VDW and bonded interactions matters.
  • before embarking on a convoluted study of a complex system, I would first do some simulations with just OPLS-AA and study the mobility of your hydroxide ions overall. I would then look into making some simulations of bulk water with KOH using ReaxFF (with manageable size systems) to see how much reactivity you will actually see on the time scale available. Based on the results of these two studies, I would decide on whether a complex study is useful. Keep in mind that the error is the worst of the two subsystems plus the error from the (mechanical) coupling. The fact that you have different parameterization balances for OPLS-AA and ReaxFF and that you have to embed a polarizable system into a non-polarizable one (which will “bleed” energy because of this inconsistency) may lead to an error that is too large to make any meaningful analysis. So it is a very high risk project and thus it is very advisable to make an assessment before jumping right in.

In addition to @akohlmey 's thoughts, you should look through the literature to see what has previously been done (being aware that many terrible MD studies can be published through peer review). This paper: https://pubs.acs.org/doi/10.1021/acs.iecr.3c00636 may give you some indication of what is feasible for MD simulations.

It is also not a good idea, especially when new, to try drawing conclusions from simulating a single system. You should start with a clear property trend from experiments, and clear competing structural hypotheses to test with your MD. (If there is only one possible structural explanation for your

For example, suppose I vary the concentration of acetate pH buffer in an ion exchange resin experiment and find that higher concentrations of acetate make the ion exchange slower. And let’s suppose either the acetate is clogging up the polymer, or acetate associates more strongly with the cations than hydroxide. Well, those are two hypotheses we can test by simulating different concentrations of acetate and hydroxide and seeing what happens to the various g(r) and diffusivities. And then even if I use a very simple molecular model, if it replicates the experimental results, then that’s good enough, especially to back up solid experimentalists.

Hi Axel,

I appreciate your suggestions greatly. I am using HPC but I think it’s still not worth consuming large resources to deal with such a system. People tend to design a small and representative system to look at the reactive stuff. The large size of my current system is more suitable for doing geometric analysis rather than probing a sophisticated ion conduction mechanism unless there exists a ‘more is different’ effect.

As I rethink my objective, maybe classical MD is enough to describe the interaction events. The idea of doing a reactive simulation seems natural but overly naive when it comes to modelling on hydroxide, not to mention in a large and complex system.

Thank you again for your suggestions.

Best,

Eiji

Hi srtee,

Thank you very much for your suggestions and the literature. I do realise there must be a matchup between the simulation and experiments and this is what I have been doing. Regarding my case, as I re-evaluate the objective, I would say the classical MD is sufficient to deal with my system. There may be no need to explicitly model the formation and breaking of hydrogen bonding

What I am asked to do now, like most other computationalists on material modelling, is to do the modelling complementary to experiments. Experimentalists sometimes underestimate the modelling complexity, which is confusing to me. But now I have figured it out in this discussion.

Best,

Eiji