Transfer from GULP FF to Lammps

GULP
1

Hello to gulp users,

We have been discovering GULP and its very impressive list of options, possibilities. Bravo.

We are trying to generate a Buckingham potential for crystalline SiS2, having at our disposal only a single crystal structure of 12 atoms (see below). To us (but we are new gulp users), the fit converges well:

**** Fit completed successfully ****
Final sum of squares = 0.000000
Final gradient norm = 0.019527

and we get the optimized parameters for the force-field. However, as we transfer both the structure and the force-field into Lammps, the simulation becomes highly unstable and diverges. We were expecting to have a gentle atomic vibration around the crystal positions. The same situation remains with additional options (relax or optim) and slightly modified parameters of the FF. Any idea of what we have been doing wrong ?
Thanks a lot
Matthieu

fit conv
cell
9.54500000 5.56400000 5.55200000 90.000000 90.000000 90.00
fractional
Si core 0.00000000 0.00000000 0.75000000
Si core 0.50000000 0.50000000 0.75000000
Si core 0.00000000 0.00000000 0.25000000
Si core 0.50000000 0.50000000 0.25000000
S core 0.38180000 0.29120000 0.50000000
S core 0.38180000 0.70880000 0.00000000
S core 0.61820000 0.29120000 0.00000000
S core 0.88180000 0.79120000 0.00000000
S core 0.11820000 0.79120000 0.50000000
S core 0.11820000 0.20880000 0.00000000
S core 0.61820000 0.70880000 0.50000000
S core 0.88180000 0.20880000 0.50000000
species 2
Si core 3.200000
S core -1.600000
end
buck
Si core Si core 0.0000 0.072000 0.0000000 0.0 8.0 0 0 0
Si core S core 1500.0 0.2512000 260.0 0.0 8.0 1 1 1
S core S core 3500.0 0.4453000 1800.0 0.0 8.0 1 1 1

dump example_essai.grs

2

Hi Matthieu
The issue you’re having is because there was insufficient information in the fit. Although the forces have gone to zero, there was nothing to constrain the curvature of the potential energy surface. If you use your fitted potentials to optimise the structure and compute the mechanical and phonon properties of your system you’ll see that you have negative elastic constants / bulk moduli and lots of imaginary phonons. This means your system is unstable and so LAMMPS is quite correct when the thing collapses in a horrible mess.
The way to avoid this is to add in curvature properties (elastic constants, phonons) in the fit. You should also try to use “relax” with the fit when your parameters are good enough to optimise the structure.
I’d also note that your C6 terms look unphysically large. I’d recommend you make them smaller and/or add a taper function to smooth the boundary as your cutoffs are short for such large C6 terms. Generally it’s also best not to fit them unless you have a lot of fitting observables over a wide range of volumes.
Hope that helps,
Julian