ReaxFF bond order cutoffs vs VMD heuristic bond assignment approach

LAMMPS LAMMPS Beginners
, ,
1

Hi everyone,

I’m conducting MD simulations of the fluorination of bilayer graphene using the ReaxFF force field. My goal is to study the system at various temperatures (700K, 500K, 300K…) to observe C-F bonding as well as interlayer C-C bonding. Recently, I’ve encountered issues with post processing my results to look for these C-C bonds. The way I do this is by using the bond information dumped in a .reax file and post processing that data using my own analysis software/script as discussed in Bond order cutoff - #4 by akohlmey. Before going further, here is my code:

variable dt equal 0.5      # Timestep in fs
variable steps_per_restart equal 1000000
variable total_steps equal 5000000

units real
atom_style charge
boundary p p f  # periodic in x and y, fixed in z

# Read graphene structure (carbon only)
read_data graphene_mini.data  

# ReaxFF potential
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield.reax.FC.txt C F 

# Define groups
group graphene type 1

# Charge equilibration
fix qeq all qeq/reaxff 100 0.0 10.0 1.0e-6 reaxff

# Prevent atoms from flying out of z bounds
fix wall all wall/reflect zlo EDGE zhi EDGE

# Temperature monitoring
compute tgraphene graphene temp

# ----------------------------
# Relaxation Phase
# ----------------------------
fix boxrelax all box/relax x 0.0 y 0.0 vmax 0.001
minimize 1.0e-6 1.0e-9 1000 10000
unfix boxrelax

# Add F atoms
region above_graphene block 1 43 1 41 8 49 units box  
region below_graphene block 1 43 1 41 -49 -8 units box 
create_atoms 2 random 1440 12345 above_graphene overlap 1.5
create_atoms 2 random 1440 67890 below_graphene overlap 1.5
group gas type 2

compute tgas gas temp

velocity all create ${T} 4928459 mom yes rot yes dist gaussian

# Soft start: Langevin thermostat + NVE
fix pre_relax all nve
fix temp_relax all langevin ${T} ${T} 100.0 123456
thermo 50
timestep ${dt}
run 3000
unfix pre_relax
unfix temp_relax

# ----------------------------
# Main Simulation Phase
# ----------------------------
fix nvt_sim all nvt temp ${T} ${T} 50

# ReaxFF bond info
fix reax_bonds all reaxff/bonds 100 bonds.reax

# Restart file every 1 million timesteps
restart ${steps_per_restart} restart.*.bin

# Dump less frequently to reduce file size
dump mydmp all custom 10000 dump.lammpstrj id type x y z q
thermo 1000

run ${total_steps}

The issue I’m facing is that I’ve also used VMD throughout my work to look for bonds since in the case of C-F bonds for example, those are displayed normally using various representations (Dynamic bonds, CPK). However, in the case of these interlayer C-C bonds, I’m getting mixed signals because the bonds aren’t shown in VMD despite the surface topology being indicative of bonding (as my graphene sheets “pucker” similarly to F-Diamane, which is the desired end goal). Below is an image of what I’m trying to describe (F atoms blue and C atoms red):

Side view of the graphene sheets showing C-F bonding and potentially missing C-C bonds (Dynamic bonds view))
Side view of the graphene sheets showing C-F bonding and potentially missing C-C bonds (Dynamic bonds view))962×521 110 KB

I’d also like to clarify that I’m talking about these interlayer C-C bonds in an indefinite manner because:

  1. the objective of my work is to investigate if these bonds will occur or not, i.e. whether they’re possible in the context of MD simulations or not
  2. Due to my lack of knowledge of a suitable bond order cutoff, I’m doubtful about the validity of my manual analysis approach and the insights I gain from that, and so I thought of similarly considering VMD as I’ve done with C-F bonds, except this time my manual analysis and VMD disagree as I’ve tried arbitrary cutoffs to explore C-C bond count trends against time and I end up getting a non-zero amount of bonds, whereas VMD shows me none.

This contrast made me take a step back and evaluate what I should “trust” more, my manual analysis (which is incomplete/not functional due to my lack of a cutoff) or VMD’s heuristic approach and what it displays/the cutoff it employs? This would be a more pleasant question if I were to know a benchmark bond order cutoff value. I’ve tried to do some research or digging on this but couldn’t quite find something helpful. However, I’m not ruling out the existence of an answer or helpful resource on the internet that I could’ve potentially missed or overlooked due to my lack of expertise with LAMMPS, VMD, and ReaxFF workflows.

Lastly, while going through the forum I ran into the following post Looking for Feedback on new Feature in LAMMPS / LAMMPS-GUI which addresses the issue of bonds not being drawn at times. I didn’t understand everything in the post but I’m curious if it can be relevant to my case at hand?

— LAMMPS version: 29 Aug 2024 - Update 1 , on Windows

Thank you in advance and apologies for any oversights or misconceptions

2

No. At least not at the moment. This feature works pretty much the same as the DynamicBonds representation in VMD by using a fixed cutoff for all bonds (which is clearly incorrect for anything beyond trivial cases where you have just covalent single bonds between elements like carbon, nitrogen, oxygen. Already hydrogen is creating a problem, so there is an additional hack (in both the VMD and the LAMMPS visualization) to disallow any bonds between two hydrogen atoms.

I have been contemplating two alternate bond algorithms:

  1. use a variable cutoff based on covalent radii of atoms (same as the VMD automatic bond guess uses). this avoids some of the problems with just one cutoff, but requires a correct identification of all elements (which is a problem in VMD, especially when feeding it LAMMPS dump files without element information).
  2. implement an interface to fix bond/reaxff and then use a bond-order cutoff instead of a distance cutoff to decide which bonds to draw.

I have been hesitant to implement either since it is a) a lot of work and b) this is “only” about visualization and not about analysis and analysis needs to be done with scripts/programs during post-processing. The current method is suitable for at least 80% of the use cases as it is.

This is completely arbitrary and thus unsuitable for a proper analysis.

Please note that any calculation with fluorine is difficult in general since already the parameterization is difficult since the underlying DFT calculations tend to struggle to represent it well.

1 Like
3

This knowledge is available from learning how ReaxFF works by studying the relevant publications with sufficient care.

Please also note, that this kind of expertise at the level you will need as a beginner in the field, you cannot obtain from an online forum. This requires to have personal supervision and tutoring by a sufficiently qualified supervisor, experienced colleague, or collaborator. In the old days (when I was a grad student), it was quite common to join the group of somebody with the necessary expertise for a few weeks/months to obtain that kind of knowledge and experience (and then share it with your colleagues locally).

1 Like