Unreasonable Benzenedithiol Structure during pulling (using reactive FF)

Dear LAMMPS users,

I am using LAMMPS 16 Feb 2016 version to run a pulling test for Benzenedithiol using reaxFF (reax/c command). Attached are the input files(data.bzdithiol, in.pull), potential and control files(ffield.reax, lmp_control), resulting force curve (Rupture_Force_Au-S.png), and structure snapshots (Minimized_structure.png and During_pulling.png). The initial structure is that I have an Au atom to mimic the tip and fix the other Au atom as the surface. The benzenedithiol (without the end H atoms) is a molecular bridge between the metal atoms. Sulfur is bonding with Au and C.

The specific question I want to ask is as the following:

After minimization, the benzenedithiol structure is planar as shown in Minimized_structure.png, which is reasonable. However, as long as I start to pull from one of the Au atom ends, the sulfur atoms starts to bend strangely as shown in During_pulling.png. But the rupture force curve makes sense with a 2.5nN rupture peak and the following noise around 0nN.

My guess for the strange conformation change of benzenedithiol is that the potential file of the reaxFF is not appropriate. I obtain this force field parameters from the supporting file of the paper “Development of Interatomic ReaxFF Potentials for Au–S–C–H Systems”.

Did anyone have similar experience about the strange Sulfur-benzene-Sulfur “dihedral”? Is that reasonable? Could any expert point out if this is because of bad LAMMPS scripts or I should dig more into the reaxFF potential parameters?

Thank you and best wishes,
Zhi

b.psf (1.4 KB)

data.bzdithiol (1.14 KB)

During_pulling.png

f.dat (372 KB)

ffield.reax.AuSCH (10.6 KB)

in.pull (3.47 KB)

lmp_control (1011 Bytes)

Minimized_structure.png

Rupture_Force_Au-S.png

There could be two reasons:

  1. Using SMD to pull just one atom is somewhat unconventional (at least to the best of my knowledge). Have you tried using velocity set or fix move to move that one Au atom?
  2. Looking at the force field file ffield.reax, Au-S-C-C dihedral (torsion) does not seem to be specifically parameterized. It uses the generic torsion interaction defined by the ‘0 1 3 0’ entry.
  3. Au-S-C angle has been parameterized, however, the equilibrium angle (theta_0) is 82.85 deg. Indicating this angle would not want to be linear.
    In sum, I would be more careful in applying this force field description for your purpose.

Ray