Change topology from GAFF to REAXff

Hi LAMMPS users,

I have simulated the mechanical and thermal physical properties of the different polymers using the GAFF force field at different temperatures below and close to the glassy state.

However, in the classical force field, it is impossible to simulate the deformation of the break of the polymer upon deformation and the formation of chemical cross-links because of the influence of temperature.

To predict breaking in deformation and formation of the crosslinking, I want to change my topology files from GAFF to REAXff. For a simple view, I need only coordinated of atoms because in REAXff there is no topology information of bonds, angles, and dihedrals. Moreover, since the chemical structure of my polymers consists of C, O, H, and N atoms, I could find the parameters in the REAXff simulation in LAMMPS to parametrize the interactions of these atoms in the model.
After performing the breaking or chemical cross-linking formation, I want to change the REAXff force field to GAFF to continue the simulation using the classical force field.

What do you think? Is there a simple way to perform such forward and backward transformations from GAFF to REAXff?

A LAMMPS data file is simply a list of particles and their properties. So you can simply delete the file’s “topology” sections (and change the header lines consistently), as well as the various “coefficient” sections. This gives you a pure set of particles’ coordinates and velocities, from which you can run your simulation.

Equally well, you can write a short LAMMPS script that reads in your previous data file (under the previous force field), runs the delete_bonds command to erase the topology (check the manual for the appropriate keywords), and write_data to accomplish a very similar result.

Either way, you can certainly run a ReaxFF simulation from the final coordinates of your GAFF simulation. Whether that is scientifically useful is something you will have to find out.

Dear @srtee

Thank you for your prompt response. So, as I understand, it is possible to convert a topology file from GAFF to ReaxFF. Backward, if I want to change the topology from ReaxFF to GAFF, I need to write three LAMMPS input files:

  1. simulation in GAFF,
  2. delete bonds, angles, and dihedrals; write ReaxFF parameters; perform minimization; and perform ReaxFF simulation (breaking bonds or formation cross-linking). Compare the results of the ReaxFF simulation with those of the experiment.
  3. write bonds, angles, and dihedrals in the GAFF force field for coordinates of the systems after simulation in ReaxFF and perform classical atomistic simulation in GAFF.

How do you think, will be possible to minimize the system after the backward transformation from ReaxFF to GAFF?

The backwards conversion can be quite difficult, since you cannot just recover the data you have used before after your system has undergone chemical reactions. Please see the discussion on the same topic here: Convert CHARGE data file to FULL data file

Dear @akohlmey

Thank you for your kind response. I also use the TopoTools plugin in VMD to create the initial atomistic configurations of the systems studied in LAMMPS. However, I had to perform quite a lot of manual work in addition to tsl scripting.

I believe that the main problem after the chemical reaction could be to assign the types of atoms and write bonds, angles, and dihedrals that can be used in GAFF. Also, the partial charges of the atoms after the chemical reaction should be changed significantly, and thus they should be recalculated by using quantum chemistry.

That is not automatically a good idea. You need to identify the molecules and figure out whether there are properly tested force field parameters for those compounds.

Have you looked at fix bond/react in LAMMPS?

ReaxFF parameters are not as universal as in quantum mechanics. Each parameter set is tuned for a specific set of compounds and thermodynamic conditions. Using those parameters for anything else than the intended purpose can likely lead to bad results.

In general, switching between force fields for parts of a complete study seems like a very bad idea. If you are going to use ReaxFF it should be suitable for all of your simulations or else you cannot trust your results.

Dear @akohlmey, for me, it is a very complex task. I believe that I could determine what substances will be created after breaking the bonds or cross-linking and then choosing the appropriate atom types to continue the simulation using GAFF.

Force field parameterization is a complex task and when you want to get dependable results you have to do what needs to be done, whether it is convenient or not. Science is merciless in that regard. If you are lacking experience you perhaps need to find a collaborator that has it.

After that I read this post Fix bond/react and How to cross-link amorphous polymer chains about the creation of crosslinking in the polymer melt, I will try to use this command first.

Thank you @akohlmey ,

I totally agree with you. Your help is very useful.

Dear @srtee and @akohlmey,

I appreciate your assistance.

The polymers I am referring to are heterocyclic polyamides containing the elements carbon (C), nitrogen (N), oxygen (O), hydrogen (H), and sulfur (S).

How can I select the appropriate REAXff parameters for the atom types in my molecules?

Only the parameters C and H were included in the primary publication by van Duin et al. in the Journal of Physical Chemistry A, 2001, volume 105, pages 9396-9409.

I found a website, specifically Interatomic Potentials Repository, that provides access to several potentials, such as the C-O-N-H potential. Nevertheless, within the REAXff documentation, I noticed more appropriate collections of parameters, such as CHONSSiNaFZr.ff. However, I am uncertain about where I may access and download these parameters.
Is it advisable to contact the developers of REAXff in order to acquire these parameters? Or there exist another open service where one could download the REAXff parameters?

You need to find a parameter set that was created for systems very similar to yours. ReaxFF parameter sets are much more specific then force fields and certainly much more specific as a low accuracy generalized force field like GAFF. Once you have the parameter file you must match the GAFF force field types with the ReaxFF parameter elements.

The first ReaxFF parameterizations are only applicable to gas phase combustion studies.

Those are usually included in the corresponding publications as supporting information. Also the group of Adri van Duin maintains a database of ReaxFF parameter sets that they know about. You need to register on their website to gain access.

Dear @akohlmey
Thank you for your help